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The  SOCRATES  contamination  model  is  described  at  length.  The  model  allows 
for  unsteady  or  steady  simulation  of  contamination  aboard  a  spacecraft  via  the 
direct  simulation  Monte  Carlo  method.  The  basis  for  the  model  is  discussed,  and 
sample  calculations  are  given  for  an  RCS  engine  firing  at  320  kilometer  altitude 
for  three  different  orientations.  These  calculations  arc  compared  to  visible 
data  in  the  6300  K  region  which  is  attributed  to  the  forbidden  transition 
between  the  0(^D)  and  0(^P)  electronic  states  of  i^cmosphcrlc  atomic  oxygon. 
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1.  OVERVIEW  OF  THE  SOCRATES  CODE 


The  direct  simulation  Monte  Carlo  method,  as  pioneered  by  G.  A.  Bird,'  provides  a 
powerful  technique  for  the  simulation  of  real  gas  flows.  It  bridges  the  gap  between  continuum  and 
free  molecular  flow,  retaining  validity  in  either  extreme.  It  can  be  used  to  describe  complex 
mixtures,  including  effects  of  chemical  reactions,  heat  conduction,  viscosity  and  diffusion  for  flows 
in  three  dimensions.  This  report  describes  the  application  of  this  technique  to  the  contamination 
problem,  considering  flow  fields  created  by  the  interaction  of  a  spacecraft  with  the  atmosphere. 
The  resultant  model  has  been  named  the  SOCRATES  code,  which  is  an  acronym  for 
^acecraft/Orbiter  Contamination  Representation  Accounting  for  Transiently  Emitted  Species. 

Contamination  of  instruments  on  space  platforms  is  an  issue  of  major  concern.  The  shuttle, 
for  instance,  gives  off  matter  through  surface  outgassing,  via  various  thrusters,  and  from  flash 
evaporators.  At  altitudes  where  the  atmospheric  mean  free  path  is  comparable  to  or  less  than 
vehicle  dimensions,  the  deposition  back  onto  instruments  will  be  largely  determined  by  the  multiple 
collision  environment  surrounding  the  spacecraft.  Even  at  higher  altitudes,  this  may  be  the 
dominant  source  of  contaminants  for  some  portions  of  the  vehicle.  In  addition  to  physical 
contaminuiion  of  shuttle  surfaces,  “radiation  contamination''  is  also  a  potential  problem  as  gases 
surrounding  the  spacecraft  collide  at  high  speed  with  atmospheric  molecules.  These  energetic 
collisions  can  lead  to  vibrational,  rotational,  and  electronic  excitation  and  subsequent  radiative 
decay.  Ions  in  tlie  vehicle  environment  may  remain  there  for  some  time  duo  to  electric  Held  forces, 
and  radiative  recombination  is  another  potential  source  of  radiation  contamination.  The  situation  is 
depicted  schematically  in  Figure  1. 

Spectral  Sciences,  Inc.  (SSI)  has  develop^  this  first  module  of  a  three-dimensional 
description  of  the  flow  field  around  a  spacecraft  so  that  contamination  can  be  accurately 
characterized  and  undersuKKt.  A  comprehensive  model  of  the  contaminant  field  surrounding  the 
space  shuttle  orhiler,  for  instance,  is  crucial  to  the  design  of  experiments  which  are  to  fly  ou  die 
shuttle  and  to  the  development  of  procedures  fur  minimizing  the  contamination.  Hie  code  is 
designed  in  a  highly  modularized  fashion  so  that  additional  physical  and  geometric  complexity  can 
be  added  as  deemed  necessary  without  requiring  major  rewriting  of  the  model.  Ihe  fitst  module  of 
SOCRATES  emphasizes  the  gaseous  interactions  and  ilieir  emis.siuns.  It  is  (mssihle  to  pul  the 
results  of  this  code  together  with  a  surface  unidei  to  gel  a  realistic  picture  of  how  a  vehicle  looks  in 
the  context  of  the  ga.suous  emission  it  produces. 

Work  has  also  progressed  on  adding  an  inner  solution  capability  to  SOCRATES,  so  that  the 
detailed  interaction  of  the  flow  fleld  with  the  vciiicle  will  be  accurately  described,  lliis  work  is 
described  in  this  report,  but  will  be  implemented  in  the  second  module  of  SOCRATES. 


ATMOSPHERIC  WIND 


Tho  basic  culcuiational  luchniquo  is  well  described  by  iis  uriginauir  in  Reference  i. 
However,  there  have  been  significant  extensions  the  ntethud  since  iltc  (uihllcaiion  of  Reference 
1.  The  present  purpose  Is  to  describe  how  tlte  technique  is  ln>plemcnied  in  SCXTRAlTiS;  but 
elemeiuary  concepts  and  relations  whicti  are  essential  to  a  colierem  explanation  are  included  here 
also. 


Tho  direct  simulation  Monte  Carlo  method  involves  storing  a  discrete  number  of  molecules 
(via  their  velocities,  positions,  and  other  pertinent  information)  in  a  computer.  The  .solution  region 
is  broken  up  into  a  number  of  separate  ceils,  and  the  solution  is  stepped  forward  in  time  in  a  two 
stage  process.  First,  tlte  molecules  are  advanced  along  their  trajectories  by  an  amount  appropriate 
to  their  velocity  and  a  time  increment.  At,„.  In  this  first  stage  some  molccule.s  will  leave  the 
solution  region,  and  some  will  be  introduced  as  determined  by  the  boundary  conditions  for  a 
particular  problem.  The  second  stage  is  to  simulate  collisions  in  eaci)  cell  appropriate  to  At,,,  so 
that  cettision  frequencies  are  pro])erly  simulated.  A  basic  hypoihcsi.s  of  the  method  is  that  if  the 
time  step  is  made  smalt  enough,  the  processes  of  translations  and  collisions  can  be  uncoupled  in  this 
manner. 
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Periodically,  the  solution  is  sampled  by  accumulating  statistical  sums  of  number  densities, 
velocities  and  other  basic  properties.  T^e  solution  is  run  repeatedly  until  statistical  deviations  are 
reduced  to  a  desired  limit,  and  then  physically  meaningful  output  quantities  are  computed  from  the 
statistical  sums.  The  number  of  molecules  represented  is  typically  many  thousand  at  a  time,  which 
is  vastly  fewer  than  the  number  occurring  in  virtually  all  real  flows.  Hence,  the  construction  of  a 
dynamically  similar  flow  to  be  simulated  in  the  computer  is  an  e.ssential  feature  of  the  method. 

The  logical  flt)w  of  the  solution  procedure  is  shown  in  Figure  2,  which  includes  the  steps 
described  above.  The  following  sections  de.scribe  in  detail  the  implementation  of  each  of  the  boxes 
shown  in  Figure  2  and  the  application  of  the  code  to  so.me  sample  problems 


DEFINE  THE  INITIAL  STATE  OF  THE  SIMULATION 


A  Diagram  of  fhc  Basic  Solution  Pr^Kcdurc  Utilizol  for  Steady  State 
Soluibns  in  the  SOCRA'i'ES  Cuntaminaiion  Model. 


2.  GRID  COORDINATES  AND  GRID  STRUCTURE 


As  discussed  in  Section  1,  the  Monte  Carlo  procedure  works  by  breaking  the  calculation 
region  up  into  cells.  A  solution  cell  should  be  a  region  in  which  no  properties  change  greatly,  i.e., 
the  dimensions  of  a  cell  should  ideally  be  small  compared  to  the  local  scale  length  of  the  flow  field. 
Collisions  are  simulated  on  a  cell-by-cell  basis,  and  molecules  can  experience  collisions  only  with 
other  molecules  in  the  same  cell.  There  is  no  other  spatial  criterion  used  for  determining  collision 
partners,  so  the  cell  determines  the  collision  environment  for  any  molecule  within  it.  In  addition  to 
defining  the  collision  environment,  the  other  major  function  of  the  cel!  structure  is  to  determine  the 
points  at  which  output  is  generated.  There  is  no  requirement  that  the  cells  be  divided  up  in  the 
same  coordinate  system  used  in  the  molecular  state  vector. 

For  Monte  Carlo  calculations,  as  for  other  types  of  computational  fluid  mechanic  analyses, 
the  selection  of  grid  geometry  is  a  critical  requirement  which  is  often  more  of  an  art  than  a  science. 
Considerations  in  the  selection  of  a  grid  are: 

•  The  grid  should  be  as  simple  as  possible,  since  the  program  must  repeatedly 
decide  which  cell  molecules  reside  in  as  they  move  throughout  the  solution 
region.  If  this  determination  required  the  solution  of  a  complex  equation  or 
sifting  through  tables,  the  entire  program  would  run  significantly  slower  than  if 
the  cell  can  be  determined  easily. 

•  The  grid  should  concentrate  cells  where  gradients  are  the  largest,  so  that  the 
least  number  of  total  cells  (and  molecules)  are  needed  to  obtain  an  accurate 
solution. 

•  The  grid  should  provide  flow  field  iiformation  where  it  is  required,  with  the 
resolution  that  is  desired  for  the  answer  of  Interest. 

The  SOCRATES  grid  structure  is  a  simple  extension  of  the  basic  Cartesian  coordinate 
system  that  is  used  elsewhere  in  the  code.  The  cells  are  determined  by  the  intersection  of  three 
families  of  planes,  each  family  being  perpendicular  to  one  of  the  coordinate  directions.  For  each 
coordinate,  there  is  a  plane  at  zero  and  subsequent  planes  proceed  outward  in  the  plus  and  minus 
coordinate  direction.  For  instance,  the  intersection  points,  Xj,  on  the  positive  x  axis  are  given  by 

•  -  :  ^  exp(B)-  i  ’ 


where  x^^^  is  the  position  of  the  last  plane  (the  edge  of  the  solution  region  in  that  direction),  N  is 
the  number  of  planes  intersecting  the  positive  x  axis,  and  B  is  an  adjustable  parameter.  For  B 
approaching  zero,  successive  planes  have  equal  spatial  increments;  and  as  B  is  increased  the  planes 
become  more  concentrated  near  the  origin.  The  same  relation  is  applied  for  planes  intersecting  the 
minus  x  direction,  with  x^^  being  replaced  by  x,„jn.  The  other  four  directions  (±y  and  ±z)  are 
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handled  in  an  analogous  fashion.  Note  that  the  B  and  N  parameters  are  specified  separately  for 
each  of  the  six  directions  away  from  the  origin,  depending  on  the  physics  of  the  problem  under 
consideration.  These  values  can  be  input  by  the  user  or  automatically  selected  by  the  program. 
The  SOCRATES  grid  structure  fulfills  the  objectives  enumerated  above  to  a  substantial  degree. 
The  relations  for  cell  boundary  locations  are  easily  inverted  to  obtain  the  cell  number 
corresponding  to  a  given  position,  and  the  parameter  of  the  distribution  allows  for  concentration  of 
cells  in  the  inner  region  while  allowing  larger  cells  further  out  where  the  gradients  are  less  severe. 

A  sample  cell  structure  resulting  from  this  technique  is  illustrated  in  Figure  3.  For  visual 
clarity,  the  number  of  planes  has  been  limited  to  two  in  each  of  the  six  directions,  since  this  is  the 
fewest  number  which  illustrates  the  uneven  spacing.  A  typical  calculation  would  have  several 
times  that  many  planes,  but  the  figure  is  difficult  enough  to  interpret  as  it  is.  A  vehicle  can  be 
arbitrarily  located  within  this  cell  structure;  though  it  should  be  located  near  the  center  of  the  grid 
structure  for  it  to  make  sense.  Similarly,  the  wind  direction  as  seen  from  the  shuttle  can  come 
from  any  direction  whatsoever;  there  is  nothing  in  the  cell  structure  or  coordinate  definition  which 
restricts  it. 


Figure  3.  A  Schematic  Showing  the  Baiic  De.sign  v>f  the  SOCRATES  Cell  Structure 
Which  Uses  Cartesian  Courdinutes  with  Uneven  Spacing. 


3.  GAS  MODEL  AND  EQUILIBRIUM  REFERENCE  PROPERTIES 


3.1  Preliminary  Equilibrium  Gas  Relations 

The  far  field  equilibrium  state  has  properties  which  are  of  relevance  to  the  flow  field 
interaction  problem  to  be  solved.  In  addition  to  serving  to  define  the  proper  outer  boundary 
conditions,  the  far  field  serves  to  define  length  and  velocity  scales  for  the  problem  which  are  then 
used  to  non-dimensionalize  the  internal  code  variables.  Even  if  this  were  not  done,  it  would 
provide  an  important  comparison  case  for  densities,  velocities,  collision  frequencies,  etc. 

For  a  rest  gas  in  equilibrium,  the  normalized  distribution  function  for  the  relative  speed,  c,, 
between  molecules  of  species  i  and  molecules  of  species  j  is  given  by^ 

fij(c,)  =  (4cja?j'^/v7)  exp(-aijc;)  .  (2) 


where 


2RoT« 


» 


(3) 


and  My  is  the  reduced  mass  of  the  pair;  i.e.. 


B 


mjmj 

nii+mj 


(4) 


with  nil  and  ntj  representing  the  ma.s$es  of  the  two  species.  In  these  relations,  T,.  is  the  far  field 
temperature  and  Rq  i.*i  the  universal  gas  constant.  (Rq  is  used  imituad  of  Boltzmann’s  constant  since 
the  molecular  masses  will  be  consistently  represented  in  atomic  mass  units  rather  than  grams.)  Ttie 
available  translational  collisional  energy  between  the  two  molecules,  E^.,  is  given  by 


(5) 


3.2  Analytical  Fonrn  of  the  Collision  Cross  Section 

Whenever  the  direct  simulation  Monte  Carlo  method  is  applied,  it  is  necos.sary  to  make 
trade-offs  between  accuracy  and  simplicity  in  molecular  miKlels.  It  dues  no  good  to  use  a  complex 
molecular  potential  surface  and  then  find  that  reasonable  computer  run  times  result  in  very  large 
stuistical  iluciuations  in  the  output.  Since  the  final  output  will  refiect  errors  in  the  statistics  as  well 
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as  errors  in  the  models,  there  is  a  strong  impetus  to  use  models  which  contain  the  essential  physics, 
but  v'hich  can  be  applied  in  a  computationally  efficient  manner.  I'he  current  state-of-the-art  is  the 
Variable-Hard-Sphsre  (VHS)  model. ^  In  this  model  molecules  have  a  collision  cross  section  which 
varies  as  an  inverse  power  of  the  available  collision  energy.  Hence,  if  Oij  is  the  collision  cross 
section  for  colUsions  of  species  i  with  species  j,  then  ffjj  is  given  by  a  relation  of  the  form 

<ry  -  A^jE^  .  (6) 

where  Ay  is  a  constant  coefficient.  It  follows  that  the  effective  diameter  for  molecules  of  species  i, 
di,  is  implicitly  defined  as  a  function  of  available  collision  energy  by  the  relation 

aa  =  Td?  =  AaE;"  .  (7) 

Ay  can  be  determined  from  a  reference  cross  section  and  velocity  via 

Au  =  •  (8) 

If  a  reference  cross  section  is  <jiven  for  a  reference  temperature  rather  than  a  reference  velocity, 
then  the  usual  choice  for  the  reference  veiocit  is  that  velocity  which  has  a  collision  energy  equal 
to  the  mean  collision  energy  occurring  in  collisions  at  the  reference  temperature.  Mathematically, 
this  is  equivalent  to 


<^r\> 

<c^a> 


(9) 


where  the  angle  brackets  indicate  averages  taken  over  ti.e  distribution  function  given  in 
Equation  (2)  evaluated  for  m^amj  and  To,=T„f.  Equation  (9)  can  be  simplified  to  give 


4(2°a>)RoT^f 

mi 


(10) 


For  simulations  involving  a  targe  number  of  .species,  t  vference  cross  sections  are  frequently 
not  available  for  all  possible,  collision  pairs.  In  ihir  case  it  is  possible  to  specify  Ay  for 
self'Collisions  only,  and  then  use  Equation  (7)  to  get  a  molecular  diameter  as  a  function  of  collision 
energy.  Tlten,  applying  the  relation 


ay  »  »l{di  +  dj)/2l2  , 


(11) 
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the  coefficient  in  Equation  (6)  for  interspecies  collisions  is  given  by 


(12) 


For  the  internal  workings  of  a  Monte  Carlo  code,  it  is  usually  more  convenient  to  express 
the  collision  cross  section  as  a  function  of  the  relative  collision  velocity  rather  than  the  collision 
energy.  This  is  simply  achieved  via  the  relation 


■2u 


(13) 


where 


By  =  Ay(Mij/2)-«  .  (14) 

The  parameter  u  can  be  related  to  ii,  the  exponent  of  distance  in  an  inverse  power 
intermolecular  force  law  via  the  relation^ 


a  = 


(15) 


Hence,  hard  sphere  molecules  (for  which  ti  goes  to  infinity)  are  represented  by  u  equal  to  zero. 
There  is  a  substantial  body  of  evidence,  however,  that  the  effective  size  of  molecules  does  indeed 
decrease  with  increasing  collision  energy,  so  a  positive  value  of  u  is  usually  a  better  choice,  a  can 
be  determined  from  molecular  beam  data,  or  from  its  macroscopic  implications.  For  example,  if  s 
is  the  exponent  for  the  variation  of  the  viscosity  coefficient  with  temperature,  then  it  can  be  shown^ 
that 

s  =  w  +  0.5  ,  (16) 

so  a  measurement  of  the  temperature  dependence  of  the  viscosity  coefficient  serves  as  an  indirect 
determination  of  u. 

In  order  to  incorporate  the  model  for  internal  energy  transfer  to  be  discujsed  in  Section  8,  it 
is  necessary  that  u  be  assumed  the  same  for  all  interactions.  This  represents  one  of  the  major 
restrictions  in  the  current  state  of  modeling. 

Although  the  sizes  of  molecules  are  allowed  to  vary  in  the  VHS  model  in  deciding  whether 
or  not  a  collision  is  to  occur,  when  a  collision  does  occur  the  post  collision  velocity  components 
are  computed  as  if  it  were  a  hard  sphere  collision  (see  Section  8).  This  results  in  a  substantia) 
computational  simplification  and  yet  retains  good  agreement  with  the  macroscopic  predictions  of 


-9- 


the  more  exact  model. ^  (See  Reference  1  for  a  discussion  of  molecular  scattering  for  general 
power  law  potentials.) 


3.3  Equilibrium  Reference  Properties  for  a  Multi-Component  Gas 

One  advantage  of  the  VHS  model  is  that  the  molecules  have  a  well  defined  cross  section,  so 
it  is  possible  to  define  a  mean  free  path  without  putting  limitations  on  the  minimum  deflection 
angle  that  is  considered.  As  is  the  general  case  for  multi-component  gases,  however,  each 
component  has  its  own  mean  free  path,  and  the  overall  mean  free  path  for  the  mixture  must  be 
defined  as  a  weighted  average  of  the  mean  free  paths  of  the  individual  species.  The  somewhat 
cumbersome  relations  required  to  calculate  the  overall  mean  free  path  are  given  here.  It  should  be 
noted  that  the  mean  free  path  is  calculated  only  once  for  a  given  problem,  so  the  computational 
effort  required  to  evaluate  it  is  completely  negligible. 

An  individual  molecule  of  species  i  will  suffer  collisions  with  molecules  of  species  j  with  a 
frequency  Py  given  by 


=  njo.<ayC,>  .  (17) 

where  is  the  number  density  of  species  j  and  <OyCr>  is  the  average  product  of  cross  section 
times  relative  velocity  for  the  two  species,  obtained  by  integrating  over  the  distribution  fiinction 
given  in  Equation  (2).  When  this  operation  is  performed,  the  result  is 


vy  =  2Bynjo,r(2  -  w)aU-‘'2/\/7 


(18) 


where  V  denotes  the  gamma  function. 

The  total  collision  frequency  for  an  individual  molecule  of  species  i,  vj,  is  obtained  by 
sununing  Equation  (17)  over  all  species,  i.e., 


J-i 


and  the  mean  free  path,  X^,  for  molecules  of  species  1  is  given  by 

Xi  “  =  “^V8RoT./(»mi)  ,  (20) 
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where  <Ci>  is  the  mean  molecular  speed  for  species  i  molecules.  The  mean  free  path  for  the  gas 
mixture,  Xa»,  is  then  defined  as  the  number  density  weighted  average  of  the  Xj  via 


^  IjS  nia»\ 


where  no#  is  the  total  number  density: 


n  - 

i=l 


(21) 


(22) 


A  useful  velocity  scale  is  given  by  v,,  defined  by 
V,  =  V2RoT«/<ra>  . 

where  <m>  is  the  reference  mean  molecular  weight,  i.e., 


(23) 


(24) 


V,  is  the  most  probable  molecular  speed  tor  molecules  of  the  mean  molecular  weight  at  the 
reference  temperature. 

3.4  Internal  Energy  Model 

The  current  state  of  modeling  for  internal  energy  effects  in  Monte  Carlo  flow  field 
simulations  is  the  phenomenological  model  of  Borgnakkc  and  Larsen.^  In  this  model,  transfer  of 
energy  between  internal  and  translational  modes  is  allowed,  but  it  is  necessary  to  assume  that  each 
species  has  a  fixed  number  of  internal  degrees  of  freedom,  This  is  equivalent  to  assuming  a 
constant  specific  heat,  Cpj,  for  each  species  which  can  be  related  to  the  number  of  internal  degrees 
of  freedom  via 


Ro 


(25) 


A',ernatively,  can  be  related  to  the  ratio  of  specific  heats  fur  species  i,  by  the  relation 
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(26) 


A  = 


5-37i 

7i-l 


The  interchange  of  internal  and  translational  energy  will  be  discussed  in  Section  8,  and  the 
selection  of  initial  conditions  will  be  discussed  in  Section  S. 
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4.  INTERNAL  REPRESENTATION 


4. 1  State  Vector 


Each  simulated  molecule  in  the  SOCRATES  code  is  represented  by  a  state  vector  which 
comprises  all  of  the  information  the  code  has  with  regard  to  that  particular  molecule.  The  state 
vector  has: 


•  Position  eiements  defining  the  location  of  the  molecule  in  Cartesian  coordinates. 

•  Three  velocity  elements,  giving  the  corresponding  velocity  components  in  the 
same  coordinate  system. 

•  A  value  for  the  internal  energy  (usually  rotational)  of  the  molecule.  Note  that 
the  basic  model  does  not  discriminate  between  internal  modes  for  a  particular 
species.  This  can  be  done,  if  desiri^,  by  introducing  separate  species  for  the 
distinct  modes. 

•  An  indicator  identifying  the  molecular  species.  This  indicator  in  turn  implies  aii 
of  the  properties  associated  with  that  species  (molecular  weight,  number  of 
internal  degrees  of  freedom,  name.  etc.). 

•  An  Indicator  giving  the  computation  cell  in  which  the  molecule  currently  resides. 
(This  could  be  c^culated  from  its  position,  but  it  is  needed  so  often  in  the 
calculation  that  the  extra  storage  location  is  justified  by  the  increase  in 
efficiency.) 

•  A  time  element,  giving  the  time  at  which  the  molecule  will  strike  a  solution 
surface  element  if  it  continues  on  its  current  irttjectory. 


4.2  Reduction  to  a  Reasonable  Number  of  Simulated  Molecules 

It  is  clearly  impossible  to  run  a  computer  simulation  with  anywhere  near  the  same  number 
of  molecules  that  exist  in  the  actual  flow  problem.  The  adjustment  that  is  made  to  make  the 
simulation  possible  is  to  artificially  increase  the  cross  section,  and  decrease  the  number  density,  by 
the  same  large  factor.  It  is  the  product  of  number  density  and  cross  section  which  determines  the 
collision  frequency  for  a  given  molecule,  and  it  is  the  collision  frequency  which  must  be  correctly 
simulated  if  the  correspondence  between  the  real  and  simulated  flows  is  to  be  accurate.  This  is  an 
essential  feature  of  the  direct  simulation  method  which  has  not  always  been  adequately  emphasized. 
It  means  tliat  the  internal  scaling  factors  do  not  proceed  on  a  strictly  dimensional  basis.  For 
example,  the  scaling  factor  for  cross  sections  is  not  the  square  of  the  scaling  factor  for  lengths. 
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4.3  Internal  Scales 


Many  problems  are  more  reasonably  handled  if  the  internal  calculations  are  carried  out  with 
scaled  or  diniensionless  values.  This  avoids  possible  problems  such  as  numerical  overflow  which 
can  cause  an  execution  time  error.  Such  errors  can  be  particularly  insidious  and  difficult  to  locate 
in  a  code  whose  very  essence  involves  the  random  combination  of  numbers. 

The  output  is  produced  in  physically  meaningful  dimensional  form.  Hence,  the  scaling  that 
is  discussed  here  is  irrelevant  (or  nearly  so)  to  the  interpretation  of  code  output;  it  is  strictly  a 
matter  of  the  internal  representation. 

The  choices  for  length  and  velocity  scales  are  Xo,  and  v,  as  defined  in  Section  3,  which  are 
used  to  non-dimerisionalize  the  position  and  velocity  elements  of  the  state  vector.  There  is  no  need 
to  provide  further  non-dimensionalization  of  mass  beyond  representing  them  in  atomic  mass  units, 
so  none  is  provided.  Hence,  the  scaling  factor  for  energy  is  Just  v;,  which  is  used  to 
non-dimensionalize  the  internal  energy  element  of  the  state  vector. 

Number  de.nsities  are  scaled  with  respect  to  the  far  field  ambient  number  density,  ng,,  which 
leaves  only  the  cross  section  scaling  factor  to  be  determined.  This  factor  follows  from  the 
condition  of  flow  similarity,  which  requires  that  the  probability  of  a  molecule  suffering  a  collision 
in  traveling  a  given  path  iength  be  accurately  simulated.  This  dimensionless  probability  can  be 
expressed  as  the  product  of  a  cross  section  times  a  number  density  times  a  path  length  (at  least  for 
small  enough  path  lengths),  and  it  is  required  that  this  product  be  the  same  for  dimensional  and 
scaled  representations.  This  implies  that  the  product  of  the  scaling  factors  for  these  three  quantities 
be  unity  and,  therefore,  that  the  cross  section  scaling  factor  be  l/(ng,Xo,).  The  internal  scaling 
factors  used  in  SOCRATES  are  summarized  in  Table  1. 


Table  1.  Scaling  Factors  Used  for  the  Internal  Representation  of 
Quantities  in  the  SOCRATES  Code.  All  Variables  are 
Defined  in  Section  3. 


PROPERTY 

SCALING  FACTOR 

Length . 

.  Xc 

Velocity  . . 

.  V, 

Time . 

.  Xo,/Vg 

Number  Density . 

Mass . 

.  a.m.u. 

Energy . 

Cross  Section . 

.  l/(n»Xg,) 
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4.4  Weighting  Factors 


*  Statistical  weighting  factors  are  a  crucial  element  of  a  successful  Monte  Carlo  simulation, 
allowing  trace  species  to  be  described  with  reasonable  accuracy.  The  weighting  factor  is  the 
number  of  "real"  molecules  that  correspond  to  each  "simulated"  molecule.  A  "simulated"  molecule 

•  corresponds  to  one  molecule’s  worth  of  storage  (one  state  vector)  allocated  in  the  program,  and  the 
weighting  factor  is  its  statistical  weight.  So,  for  example,  the  total  number  density  in  a  cell,  nggn 
can  be  expressed 

Hceii  =  2^—  ’ 
i=I 

where  Nj  indicates  the  number  of  simulated  molecules  of  species  i  in  the  cell,  W;  is  the  weighting 
factor  for  that  species  in  that  cell,  V  is  the  cell  volume,  and  p  is  the  number  of  species.  The 
product  NjWj  that  appears  in  Equation  (27)  is  termed  the  number  of  "real"  molecules  of  species  i  in 
the  cell.  Note  that  n^QU  as  calculated  by  Equation  (27)  is  a  scaled  value;  it  would  have  to  be 
multiplied  by  ng,.  as  shown  in  Table  I,  to  become  a  dimensional  evaluation  of  the  number  density. 

The  weighting  factors  used  in  SOCRATES  are  dependent  on  cell  and  species.  Hence,  flow 
fields  where  a  given  species  is  much  more  dominant  in  one  portion  of  the  solution  region  than 
another  can  be  accurately  represented. 

A  critical  error  that  can  occur  in  Monte  Carlo  codes  is  to  have  the  number  of  simulated 
molecules  exceed  the  dimensioned  limit  of  the  code.  On  the  other  hand,  it  is  generally  desirable  to 
have  as  many  molecules  as  is  feasible  to  obtain  good  statistics.  Resolution  of  these  conflicting 
considerations  is  compiicated  by  lack  of  a  priori  knowledge  of  what  the  species  number  densities 
will  be  as  a  function  of  space  and  time.  The  way  the  resolution  is  achieved  is  by  a  dynamic 
adjustment  of  the  weighting  factors,  as  required.  This  keeps  the  number  of  simulated  molecules 
more  or  less  constant  while  allowing  the  number  of  real  molecules  to  adjust  as  the  solution  evolves. 
The  introduction  of  weighting  factors,  with  the  ability  to  adjust  them  us  the  solution  demands,  is  an 
Important  feature  of  a  Monte  Carlo  simulation  which  is  to  be  usable  by  non-experts. 
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5.  INITIAL  AND  BOUNDARY  CONDITIONS 


5.1  Initial  Conditions 

Since  the  direct  simulation  Monte  Carlo  method  is  inherently  an  unsteady  technique,  an 
initial  state  must  be  specified  in  order  to  advance  the  solution.  For  situations  where  a  steady  state 
result  is  desired,  it  is  obtained  as  the  long-time  solution  to  an  unsteady  problem.  In  this  case  the 
initial  conditions  have  no  effect  on  the  eventual  solution,  but  they  may  well  have  an  impact  on  the 
speed  with  which  that  state  is  achieved.  For  steady  state  solutions,  SOCRATES  simply  starts  with 
an  evacuated  solution  region.  For  unsteady  solutions,  however,  it  is  necessary  to  start  with  a 
molecular  distribution  which  is  representative  of  the  conditions  at  the  start  of  the  desired 
simulation.  For  SOCRATES  these  conditions  correspond  to  a  uniform  flow  with  the  translational 
and  internal  modes  being  in  equiiibrium.  The  specification  of  the  initial  conditions  for  unsteady 
runs,  therefore,  involves  determining  the  state  vector  elements  consistent  with  this  condition  for  the 
desired  number  of  molecules. 


5.1.1  Number  of  Simulated  Molecules  and  Their  Weightini;  Factors 

The  desired  number  of  simulated  molecules  fur  each  species  in  each  cell,  M^,  is  an  input 
quantity,  (Typically,  simulations  aim  for  a  total  number  of  molecules  per  cell  in  the  neighborhood 
of  20.)  Given  the  Initial  number  density  to  be  simulated  for  a  species,  nj,  (which  will  have  been 
automatically  converted  to  internal  dimensions  -  see  Section  4)  the  weighting  factor  for  qiecies  i  in 
a  given  cell  is  simply 

Vn, 

'*'1  =  w  • 


where  V  is  the  cell  volume.  If  a  species  is  nut  initially  present  in  a  cell,  then  it  is  assigned  an 
initial  weighting  factor  of  zero.  If  simulated  molecules  come  into  the  cell,  the  weighting  factor 
from  their  place  of  origin  will  be  used  to  initialize  the  weighting  factor  in  the  cell.  As  the  .solution 
proceeds,  the  weighting  factors  are  automatically  adjusted  to  keep  the  average  number  of  simulated 
molecules  of  each  species  in  each  ceil  approximately  equal  to  M«.  The  one  exception  to  this  rule  is 
that  the  code  will  always  try  to  keep  a  sufficient  number  of  major  species  in  a  cell  to  guarantee  a 
prqier  collision  environment,  and  this  sometimes  requires  more  than  M,.  major  species  molecules. 


5.1.2  Initial  Po.sitions 

The  initial  molecules  assigned  to  a  cell  should  have  an  equal  probability  of  being  placed  in 
any  volume  clement  of  tiie  cell.  For  the  hexahedral  celts  of  SOCRATES,  this  simply  involves 
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selecting  each  of  the  position  elements  at  random  from  the  range  appropriate  to  the  cell  in  question. 
That  is,  the  x  position  is  selected  via  the  equation 


X 


+  ^(X, 


mux  ’ 


(29) 


where  x^i„  and  x^  are  the  positions  of  the  x-faces  of  the  cell  in  question,  and  0  denotes  a  random 
variable  which  has  equal  probability  of  lying  anywhere  in  the  interval  {0,1}.  (/J  will  appear  often 
in  this  report.  Each  appearance  corresponds,  of  course,  to  a  unique  evaluation  of  the  random 
variable.)  The  other  position  elements  are  selected  analogously. 


5.1.3  Initial  Velocity  Components 

The  thermal  velocity  components  for  a  molecule  in  translational  equilibrium  (neglecting,  for 
the  moment,  any  mean  ftow  contribution)  should  be  selected  from  a  normalized  Maxwellian 
velocity  distribution,  fQ(v),  given  by 

fo  =  a  exp|-(ov)-l/‘/ir  ,  (30) 


where 


a  «  Vm/(2RoTo.)  ,  (31) 

m  is  the  species  molecular  weight,  Rq  is  the  universal  gas  constant  and  T„  is  the  temperature. 
Equation  (31)  applies  for  each  of  the  molecular  velocity  components  and  mu.st  be  sampled  three 
times  for  each  molecule  that  comprises  the  initial  state  of  the  simulation.  A  method  for  directly 


sampling  from  this  distribution  is 

A,  =  , 

(32) 

Aj  ®  “V'-iogO?)  . 

(33) 

V  «  AisinCAj)  . 

(34) 

After  the  thermal  velocity  components  are  determined  for  each  molecule,  then  any  mean  ftow 
velocity  is  simply  added  on.  The  velocities  are  then  tran.vformed  to  internal  units. 
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5.1.4  Initial  Internal  Energies 


The  only  remaining  element  of  the  state  vector  to  be  specified  is  the  internal  energy. 
Internal  energies  for  a  gas  in  equilibrium  are  distributed  according  to  the  normalized  distribution 
function  fj  given  by 


‘  2f'2r(fy2) 


(35) 


where  f  represents  the  number  of  internal  degrees  of  freedom  for  the  species  in  question.  F  is  the 
gamma  function,  and  (  is  a  dimensionless  internal  energy,  i.e.. 


^  RoTo. 


(36) 


where  E|  is  the  internal  energy.  Equation  (35)  is  a  representation  of  the  chi  square  distribution  for 
f  degrees  of  freedom;  procedures  for  santpling  from  this  distribution  are  given  in  Reference  5. 


5.2  Source  Boundary  Conditions 

The  introduction  of  source  molecules  into  the  simulation  is  a  tmundary  condition  which 
depends  on  the  specific  model  for  the  source  in  question.  SOCKAI  ES  includes  a  *core  tlow* 
source,  which  describes  the  contamination  that  results  from  the  scattering  of  the  tlow  from  a 
thruster  back  onto  the  shuttle.  For  this  source,  it  is  important  to  have  a  description  of  the  main 
exhaust  flow  away  fr<Mn  a  thruster.  (This  is  to  be  contrasted  with  a  source  which  describes  the 
direct  contamination  of  surfaces  via  impact  of  exhaust  gases.  Since  the  thrusters  are  not  poini«xl 
directly  at  shuttle  surfaces,  it  is  the  small  portion  of  the  flow  which  leaves  the  thruster  m  a  large 
angle  from  the  exhaust  cemeriine  which  is  important  in  this  case.) 

The  plume  gases  expand  upon  leaving  the  exit  plane  and  adopt  an  essemiidly  radial  flow 
profile  over  a  distance  which  is  on  the  order  of  exit  plane  dimensiotu.  Since  this  distance  is  small 
compared  to  the  length  scales  of  the  interaction  of  the  plume  with  the  atmosphere,  it  is  appropriate 
to  replace  the  noaate  by  a  point  source  of  exhaust  molecules  traveling  at  their  thermodynamic 
timlUng  speed  with  an  uodisturbi^  number  density  distribution  given  by 

he  »  ,  (37) 

whiere  D  gives  the  axial  number  density  decay  and  0  represents  the  angle  from  the  thrust  axis.  The 
r  appearing  in  Equation  (37)  is  the  spherical  radius,  giving  the  total  distaiwc  from  the  source.  Ihe 
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particular  form  for  f(5)  that  is  used  in  SOCRATES  is  an  asymptotic  form  of  that  proposed  by 
Brook,*  namely 

f(^)  =  exp{-X“(I -cos(e)]}  ,  (38) 

where 


(39) 

(40) 

(41) 

(42) 


In  the  above  relations,  u^,  n^  and  denote  the  exit  velocity,  Mach  number,  number 
density  and  area,  respectively;  is  the  nozzle  divergence  half  angle,  and  u,„  and  7  are  the 
thermodynamic  limiting  speed  and  the  ratio  of  specific  heats. 

In  a  solution  time  step  At^,  njU,AeAi,„  real  molecules  are  introduced.  Each  molecule  is 
assigned  an  angle  $  which  is  chosen  to  be  consistent  with  Equation  (38)  via 

e  =  cos-‘(l  +  C,iog(l  -q/3)|  ,  (43) 

where 


and 

C2  »  1  -  exp(-2X2) 


(44) 


(45) 
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An  azimuthal  angle  is  selected  at  random  for  the  molecule,  and  then  the  resulting  velocity  is 
represented  in  the  basic  Cartesian  coordinates  utilized  by  SOCRATES.  The  molecuie  is  then 
advanced  from  the  source  appropriate  to  a  speed  of  u,n  and  a  time  increment  which  is  a  random 
fraction  of  At^.  The  process  is  repeated  until  the  proper  number  of  simulated  molecules  of  each 
exhaust  species  have  been  introduced. 


5.3  Atmospheric  Boundary  Condition 

The  atmospheric  boundary  condition  for  SOCRATES  is  that  molecules  should  be  introduced 
into  the  solution  region  from  the  outer  boundaries  in  such  a  way  as  to  simulate  the  undisturbed 
ambient  flow  outside  of  hie  solution  region. 


5.3.1  Molecular  Flux  Across  a  Surface  Element 


The  relations  for  molecular  flux  across  an  infinitesimal  surface  element  are  given  in 
Reference  1.  If  q  is  the  molecular  flux  (molecules  per  unit  area  per  unit  time)  crossing  a  g.ven 
surface  element,  then  q  is  given  by 

tto»  ,  .  <t.  .  .  .  . 

q  «  -^{exp^-w'j/v  vr  +  wn+en(w)i)  ,  t40) 


where 


A  =  Vm/2RoT« 


and 


(47) 


W  =  AUa,COS(0)  .  (48) 

In  these  relations,  ng,  and  T,,  represent  the  ambient  number  density  and  temperature  fur  the 
species  in  question;  m  represents  its  molecular  weight;  Uq,  represents  tiie  mean  flow  velocity:  and  & 
is  the  angle  between  the  inward  surface  normal  and  the  mean  flow  direction.  Ihe  flux  given  by 
Equation  (46)  is  positive  for  all  values  of  6,  reflecting  the  distribution  of  molecular  velocities. 
However,  it  does  become  exponentially  small  fur  large  negative  w.  (Note  that  these  relations  must 
be  applied  on  a  species-by-species  basis;  each  species  has  a  different  spread  in  its  velocity 
distribution  by  virtue  of  its  different  molecular  weight.) 

The  application  of  Equation  (46)  to  the  flat  surfaces  comprising  the  SOCRATES  outer 
boundary  is  direct,  since  0  does  not  change  along  the  face.  The  total  number  of  molecules  to 
introduce  for  a  time  step  of  At,„  is  simply  qAt,„A,  where  A  is  the  area  of  the  flat  face  in  question. 
Since  the  flux  is  constant  over  a  flat  face,  each  position  on  the  face  is  equally  likely  as  a  point  for 
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molecular  entry.  Hence,  for  a  flat  face,  the  starting  molecular  position  can  be  simply  obtained  by 
selecting  a  point  at  random  on  the  face. 


5.3.2  Incoming  Molecular  Velocity  Components 

For  each  molecule  that  is  introduced,  a  local  orthogonal  coordinate  system  is  set  up  such 
that  one  direction  is  in  the  direction  of  the  inward  surface  normal.  Velocity  components  are  first 
determined  in  terms  of  this  local  coordinate  system  and  then  transformed  to  the  main  code 
coordinate  system.  In  the  local  coordinate  system,  the  velocity  components  parallel  to  the  surface 
are  determined  as  discussed  above  for  molecules  in  the  initial  condition.  The  inward  component  of 
velocity  must  be  selected  in  proportion  to  the  distribution  h(v)  given  by 

h(v)  =  V  exp{-(Q((v  -  <v>)l-}  ,  (49) 

where  <v>  is  the  component  of  the  mean  flow  velocity  in  the  inward  normal  direction  and  a  is  as 
given  in  Equation  (31).  It  is  possible  for  <v>  to  be  negative,  but  all  incoming -molecules  must 
have  a  positive  v  value  by  definition.  Hence,  this  distribution  is  only  sampled  for  positive  values 
of  V.  The  sampling  is  done  via  the  acceptance-rejection  technique. 


6.  MOLECULAR  TRANSLATIONS 


As  discussed  in  Section  1,  an  essential  element  of  the  direct  simulation  Monte  Carlo  method 
is  the  periodic  advancement  of  simulated  molecules  along  their  trajectories.  Formally,  this  is 
accomplished  by  updating  the  position  and  velocity  elements  of  the  state  vector.  The  specific 
procedures  for  doing  this  depend  on  the  coordinate  system  in  which  the  state  vector  elements  are 
represented. 

6.1  Molecular  Translations  in  Cartesian  Coordinates 

In  Cartesian  coordinates,  the  translation  is  very  direct.  Let  x,  y,  and  z  represent  the 
position  coordinates  and  u,  v,  and  w  the  corresponding  velocity  coordinates.  If  initial  and  final 
values  of  the  state  vector  are  represented  by  a  0  and  1  subscript  respectively,  then  the  updated  state 
vector  elements  corresponding  to  a  translation  through  a  time  step  At  are  given  by 


Xi  =  Xq  +  UoAt  , 

(50) 

yt  =  yo  +  voAt  , 

(51) 

Z,  =  Zq  +  WqAI  , 

(52) 

Uj  =  Uq  , 

(53) 

< 

II 

< 

o 

(54) 

u 

< 

o 

(55) 

6.2  Molecular  Cloning 

When  a  simulated  molecule  is  translated  from  one  cell  to  another,  the  weighting  factor  tor 
that  species  will  generally  be  different  in  the  new  cell.  Since  it  is  the  number  of  real  molecules 
rather  than  the  number  of  simulated  molecules  which  must  be  preserved  when  crossing  cell 
boundaries  (statistically,  at  least),  it  is  necessary  to  correct  fur  the  distinct  weighting  factors  (see 
Subsection  4.4). 

If  the  weighting  factor  before  translation  is  Wq,  then  the  simulated  molecule  represents  that 
many  real  molecules.  If  the  weighting  factor  in  the  new  cell  is  W|,  then  Wq/Wi  simulated 
molecules  would  be  required  to  represent  the  same  number  of  real  molecules  in  the  new  cell.  If 
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this  ratio  were  a  whole  integer,  then  this  could  be  accomplished  by  introducing  that  many  "clones" 
of  the  simulated  molecules  in  the  new  cell.  That  is,  Wq/Wi  simulated  molecules  would  be  placed 
in  the  new  cell,  ail  with  the  same  state  vector. 

When  the  number  WqAVj  is  not  an  integer  (the  usual  case,  of  course),  then  the  cloning  must 
be  done  on  a  statistical  basis.  So,  for  instance,  if  Wq/Wj  were  equal  to  2.7,  then  30%  of  the  time 
two  clones  would  be  produced,  and  70%  of  the  time  three  clones  would  be  produced.  Note  that 
the  ratio  may  be  less  than  unity,  and  the  molecuie  may  not  be  introduced  into  the  new  cell  at  all. 
(In  which  case  the  molecule  is  removed  from  the  simulation.) 

Cloning  is  a  necessary  evil  inherent  in  a  system  with  spatially  varying  weighting  factors.  It 
enables  such  a  system  to  maintain  the  statistically  correct  flux  of  mass  and  momentum  across  cell 
boundaries,  but  it  misrepresents  the  flux  of  randomized  or  thermal  energy.  This  can  be  seen  by  an 
extreme  case  where  a  very  large  number  of  clones  is  produced  when  a  simulated  molecule  crosses  a 
cell  boundary.  The  resulting  molecules  in  the  new  cell  have  the  correct  mass  and  momentum  flux, 
but  since  they  all  have  precisely  the  same  velocity  they  have  a  null  relative  velocity  and,  therefore, 
a  zero  temperature.  If  the  weighting  factors  are  not  too  different  between  adjacent  cells,  then  the 
errors  introduced  by  this  process  are  acceptably  small.  However,  it  does  mean  that  one  cannot 
arbitrarily  improve  .statistics  in  one  portion  of  the  solution  region  by  selectively  reducing  the 
weighting  factors  there.  This  was  a  difficulty  which  was  encountered  in  the  eariy  stages  of  the 
direct  simulation  Monte  Carlo  method  while  trying  to  improve  statistics  along  the  axis  of 
axisymmetric  simulations,  since  the  cell  volumes  (and,  therefore,  the  sample  sizes)  tend  to  be 
smallest  on  the  axis. 

As  is  also  the  case  for  simulated  molecules  produced  via  chemical  reactions,  it  is  possible 
for  the  weighting  factors  between  successive  cells  to  be  so  different  that  a  prohibitively  large 
number  of  simulated  molecules  would  be  required  to  produce  the  same  number  of  real  molecules. 
The  codes  sense  when  a  disproportionate  number  of  simulated  molecules  are  being  produced  for  a 
given  species  and  ceil  and  adjust  the  weighting  factor  automatically.  As  the  weighting  factor  is 
increased,  a  proportionate  fraction  of  molecules  of  that  species  and  cell  are  removed  from  the 
simulation  in  order  to  keep  the  number  of  real  molecules  properly  represented.  This  process 
enables  the  weighting  factors  to  seek  their  own  proper  level  without  a  priori  knowledge  of  the 
solution.  (Periodically,  the  cells  are  examined  to  determine  if  a  certain  species  has  been 
underrepresented  in  terms  of  its  number  of  simulated  molecules.  If  this  is  found  to  be  the  case, 
then  the  weighting  factor  is  decreased,  allowing  weighting  factors  to  float  downwards  as  well  as 
upwards.  It  is  the  danger  of  weighting  factors  being  too  small,  causing  an  overflow  of  code 
dimensions,  which  is  most  critical,  however.) 


% 
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7.  COLLISION  SAMPLING  IN  A  MULTI-COMPONENT  VHS  GAS 


7.1  General  Considerations  and  Approach 

The  two  general  considerations  in  the  sampling  of  collisions  are,  as  usual,  accuracy  and 
efficiency  of  the  simulation.  As  far  as  accuracy  is  concerned,  it  is  crucial  that  the  method  in  which 
molecules  are  selected  for  collisions  be  proper.  The  correct  collision  frequency  must  be  simulated 
between  various  species  and,  in  fact,  between  the  different  portions  of  the  velocity  phase  space  for 
the  various  species.  Furthermore,  this  frequency  of  simulated  collisions  must  remain  correct 
without  any  requirements  put  on  the  velocity  distribution  function;  it  certainly  must  not  be  assumed 
that  there  is  a  Maxwellian  velocity  distribution. 

As  far  as  efficiency  is  concerned,  it  is  highly  desirable  to  use  a  method  of  collision  sampling 
involving  a  computational  effort  which  is  proportional  to  the  number  of  simulated  molecules,  N,  in 
a  cell.  Methods  which  are  proportional  to  a  power  of  N  greater  than  unity  can  become 
prohibitively  time  consuming  as  the  number  of  molecules  is  increased  -  a  limit  which  should  be 
made  as  accessible  as  possible  for  obvious  physical  reasons. 


7.2  Collisiuu  Stiiupliii^  fot  a  Single  CorApouent  Gas 

The  simplified  situation  of  a  simulation  involving  only  one  species  is  considered  here.  This 
problem  is  significant  in  part  due  to  all  the  attention  it  has  received  and,  as  will  be  seen,  it  serves 
as  an  important  reference  case.  When  there  is  just  one  species,  then  there  is  just  one  gas  kinetic 
cross  section  (though  it  is  still,  of  course,  a  function  of  collision  energy),  just  one  molecular  weight 
and  just  one  weighting  factor  for  each  cell.  In  short,  just  one  of  everything  that  has  a  molecular 
subscript.  Hence,  in  this  subsection  ail  such  quantities  will  be  presented  without  subscripts.  The 
most  important  simplification  of  having  a  single  species  is  that  there  is  just  one  collision  class,  i.e., 
only  self-collisions  of  the  given  species  with  itself  are  possible. 


7.2.1  Collision  Pair  Selection 


As  discussed  in  Section  1,  collisions  are  sampled  on  a  cell-by-cell  basis  until  the  number  of 
collisions  simulated  is  appropriate  to  the  overall  solution  time  step,  At,,,.  The  only  spatial 
requirement  placed  on  potential  collision  partners  is  that  they  be  within  the  same  cell.  In 
particular,  it  is  not  required  that  they  be  within  a  molecular  diameter  of  each  other.  (Note  that  if 
all  pairs  of  molecules  were  inspected  to  find  those  that  were  sufficiently  close  to  each  other,  this 
would  involve  a  computational  effort  in  proportion  to  the  square  of  the  number  of  molecules  in  the 
cell.)  The  rationale  for  this  is  that  the  cells  should  be  small  enough  so  that  macroscopic  properties 
can  be  assumed  constant  across  the  cell.  When  this  is  the  case,  then  a  molecule  within  the  cell  can 
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be  considered  typical  of  a  moiecule  which  might  exist  anywhere  within  the  cell,  and  molecular 
location  can  be  ignored  when  selecting  potential  collision  pairs. 

Spatial  consideration  aside,  the  probability  of  any  two  molecules  experiencing  a  collision  is 
proportional  to  aCp  the  product  of  their  mutual  cross  section  times  their  relative  velocity.  This 
probability  is  correctly  simulated  via  an  application  of  the  acceptance-rejection  technique.  A 
maximum  value  for  aCp  (ffCf)^,  is  stored  for  each  cell.  (This  value  is  updated  whenever  a  greater 
value  is  encountered.)  Pairs  of  molecules  are  selected  at  random  from  the  cell,  and  ac^  for  that  pair 
is  calculated.  The  ratio  r,  defined  by 


gCr 


(56) 


is  determined.  A  random  variable,  /3,  is  then  generated,  and  the  pair  of  molecules  is  accepted  as 
collision  partners  if  r  is  greater  than  /3.  This  produces  the  proper  relative  collision  probability 
without  regard  to  the  exi.sting  velocity  distribution  function. 


7.2.2  Collision  Time  Counter  for  a  Single  Component  Gas 

The  volumetric  collision  frequency  for  a  single  component  gas,  v,  (collisions  per  unit 
volume  per  unit  time),  is  given  by 


V 


(57) 


where,  as  in  Section  3,  n  represents  the  number  density  of  the  species,  and  <orc,>  is  the  average 
product  of  collision  cross  section  and  relative  velocity.  At  first  inspection,  it  would  seem  from 
Equation  (57)  that  a  correct  simulation  of  collision  frequency  would  require  evaluation  of  <aCf>, 
which  would  mean  that  all  pairs  of  molecules  in  a  cell  would  have  to  be  considered.  Such  a 
procedure  involves  a  computational  effort  proportional  to  N'  and  is  to  be  avoided,  if  possible,  in 
preference  to  a  method  which  is  proportional  simply  to  N.  The  alternative  approach,  introduced  by 
Bird,^  is  the  time  counter  approach.  For  each  collision  a  time  counter,  t^,  is  incremented  by  an 
amount  which  depends  on  the  relative  velocity  of  the  collision.  Collision  sampling  continues  in  a 
cell  until  its  time  counter  has  been  advanced  beyond  the  overall  flow  simulation  time,  at  which  time 
the  code  proceeds  to  the  next  cell  (which  has  its  own  time  counter).  The  time  counter  increment, 
At^,  is  given  by 


(58) 
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where  V  is  the  cell  volume  and  n  is  the  species  number  density  given  by 


! 


(59) 


with  W  being  the  weighting  factor  for  the  species.  It  should  be  stressed  that  Equation  (58)  applies 
for  each  real  collision.  As  is  discussed  in  Subsections  4.4  and  8.4,  each  simulated  collision 
corresponds  to  W  real  collisions,  so  when  a  simulated  collision  occurs  the  actual  applied  increment 
to  t,.  is  W  times  the  value  given  by  Equation  (58).  A  demonstration  of  the  validity  of  Equation  (58) 
is  given  in  Reference  7. 


7.3  Collision  Class  Sampling  in  Gas  Mixtures 


The  above  procedure  for  a  single  species  gas  can  be  extended  to  a  multi-component  mixture 
via  consideration  of  distinct  collision  classes.  In  this  approach,  collision  cla.sses  are  defined  by  the 
colliding  pair  identities.  Hence,  if  there  are  p  species  in  the  simulation  then  there  are  p(p+l)/2 
collision  classes,  which  can  be  identified  by  the  subscripts  of  the  corresponding  molecular  pair. 
(The  number  of  classes  is  not  p^  since  the  order  of  molecule  specification  is  not  taken  to  matter  in 


deiermuilng  a  collision  class.  Hence,  the  class  Idcntincd  by  the  subscripts  1,J  is  not  disti.net  from 
the  class  identified  by  the  subscripts  j,i.) 


In  collision  class  sampling  each  collision  class  is  sampled  separately,  and  the  collision 
sampling  in  a  cell  is  not  complete  until  all  classes  have  been  considered.  Each  collision  class  has 
Its  own  stored  value  of  (uyCf)ny,x  and  its  own  separate  time  counter,  t^jj.  It  can  be  shown  that  the 
appropriate  time  counter  increment  in  this  case  is 


njnjVayC,  * 


(60) 


where  is  the  Kronecker  delta  which  is  unity  for  i-J  and  zero  otherwise.  As  in  the  previous 
section,  the  above  increment  applies  for  each  real  collision.  A  simulated  collision  usually 
corresponds  to  real  collisions,  where  Wl  Is  the  lesser  of  Wj  and  Wj  (see  Subsection  8.4),  so 
when  a  simulated  collision  occurs,  the  applied  increment  to  t^y  is  times  the  result  of 
Equation  (60). 
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7.4  Global  Collision  Sampling  in  a  Gas  Mixture 

Although  the  procedure  described  above  is  quite  reasonable  for,  say,  a  two-component 
mixture,  it  becomes  exceedingly  complicated  as  the  number  of  species  increases.  For  10  species, 
for  instance,  the  program  must  loop  over  55  distinct  collision  classes  for  each  cell,  and  storage 
must  be  allocated  for  110  quantities  in  each  cell.  As  the  number  of  species  increases,  the  storage 
requirement  for  the  collision  sampling  constants  quickly  becomes  greater  than  the  storage  required 
for  the  molecular  state  vectors!  The  obvious  simplification  is  to  search  for  a  technique  where 
collisions  are  simulated  simultaneously  for  ail  collision  classes,  with  each  class  having  its  proper 
relative  probability  of  being  selected.  The  overall  collision  sampling  then  continues  untii  a  single 
time  counter  indicates  that  sufficient  collisions  have  been  sampled  in  the  current  time  step  and  cell. 


7.4.1  Global  Collision  Time  Counter 


If  molecular  pairs  are  selected  for  collisions  such  that  the  various  collision  classes 
automatically  appear  with  the  proper  relative  frequency  (see  below),  then  it  is  not  necessary  to 
consider  separate  time  counters  for  all  the  various  collision  cla-sses.  One  approach  that  could  then 
be  applied  is  to  keep  a  collision  time  counter  for  just  one  collision  class  and  increment  it  when 
collisions  of  that  class  occur.  If  the  various  collision  classes  are  being  selected  according  to  their 
correct  relative  frequency,  then  simulating  the  proper  frequency  for  one  collision  class  will  ensure, 
in  the  long  run,  that  aii  coiiisiun  classes  are  occulting  with  ihe  cuiieci  tTcqueiiCy.  A  disauvantage 
with  this  approach  is  the  necessity  of  making  an  arbitrary  choice  for  the  collision  class  which  is  to 
have  a  time  counter.  Furthermore,  there  may  be  no  good  choice  for  a  reacting  flow  where  the 
dominant  species  can  vary  strongly  from  place  to  place.  (Clearly,  one  would  not  want  to  select  a 
class  of  collision  that  does  not  occur  in  a  given  cell,  since  the  result  would  be  a  never-ending 
sampling  of  collisions  of  other  classes.) 

The  preferred  approach  is  to  define  a  global  collision  time  counter,  t^,  which  is  a  weighted 
average  of  the  time  counters  of  all  collision  classes;  i.e.. 


(61) 


j»t 


where 


c  =  £  1:°. 

••>1  j“i 


(62) 


and  the  are  non-negative  coefficients  which  can  be  selected  at  will.  Note  that  in  this 
formulation  every  collision  class  will  result  in  some  increment  of  the  global  time  counter  (unless 
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Djj  is  zero  for  that  class),  so  the  collision  sampling  frequency  is  not  dependent  on  any  one  collision 
class. 


It  remains,  of  course,  to  specify  the  Dy.  A  very  convenient  choice  is  given  by 


"inj 

1  +  5ij 


(63) 


Firstly,  Equation  (63)  is  convenient  because  it  tends  to  make  the  collision  classes  with  the  higher 
collision  frequencies  count  more,  resulting  in  good  statistics  for  tg  irrespective  of  cell  location. 
(Note  that  Djj  is  cell  dependent  since  the  species  number  densities  are  cell  dependent.)  Secondly, 
Equation  (63)  results  in  a  particularly  convenient  form  for  tg.  The  normalization  factor  given  in 
Equation  (62)  can  be  summed  analytically  to  give 


(64) 


Hence,  a  collision  of  class  ij,  which  would  produce  an  increment  of  At^ij  to  its  own  time  counter 
produces  an  increme.nt  Atg  to  tg  given  by 


At, 


n2(l  +  e^) 


Atoy  , 


(65) 


where,  again,  n  is  the  total  number  density  of  all  species  in  the  cell.  If  Equation  (60)  is  substituted 
into  Equation  (65),  the  result  is 


2 

Vn2ayc, 


(66) 


Equation  (66)  is  extremely  significant  since  it  recaptures  the  precise  form  of  the  time  counter 
increment  for  a  single  species  (Equation  (58)),  but  indicates  that  it  is  completeiy  vaiid  for  a 
multi-component  mixture  so  iong  as  the  various  collision  classes  are  sampled  with  the  proper 
relative  frequency. 


7.4.2  Collision  Pair  Selection  In  Multi-Component  Mixtures 

When  considering  selection  of  collision  pairs,  it  is  cruciut  to  remember  the  distinction 
between  real  and  simulated  molecules  discussed  in  Subsection  4.4.  Given  two  simulated  molecules 
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selected  at  random  from  within  the  cell,  the  probability  of  their  having  a  real  collision  is 
proportional  to  WiWjajjCp  However,  real  collisions  cannot  happen  individually;  they  come  at 
a  time,  where  Wl  is  the  lesser  of  Wj  and  Wj.  Hence,  when  a  collision  is  decided  upon  in  the 
program,  of  them  will  occur.  To  compensate  for  this,  potential  collision  pairs  should  be 
accepted  for  a  collision  according  to  the  size  of  Q  given  by 

Q  =  Wuayc,  ,  (67) 

where  is  the  greater  of  W;  and  Wj.  The  relative  frequency  of  real  ij  collisions  will  then  be 
proportional  to  the  product  QWl  (the  relative  probability  of  a  pair  being  accepted  for  a  collision 
times  the  number  of  real  collisions  occurring  when  the  pair  is  accepted),  which  is  the  desired 
relation.  Selection  of  collision  pairs  with  the  correct  relative  frequency  then  assures  that 
incrementing  the  global  time  counter  as  discussed  above  will  give  a  statistically  correct  sampling  of 
ail  collision  classes  simultaneously. 


7.4.3  Summary  of  Collision  Sampling  in  Multi-Component  Mixtures 

The  results  of  this  subsection  can  be  summarized  via  the  following  procedure  for  the 
sampling  of  collisions; 

•  Each  ci'l!  has  a  (current)  value  of  Q,^,  the  largest  value  of  Q  which  has  been 
encountered  so  far  In  the  collision  sampling  process.  Whenever  a  larger  value 
of  Q  is  encountered,  is  set  equal  to  that  larger  value. 

•  Each  cell  has  a  current  value  of  the  global  time  counter,  tg. 

•  Pairs  of  simulated  molecules  are  selected  at  random  from  all  molecules  within 
the  cell. 

•  For  each  pair.  Q,  (as  defined  by  Equation  (67))  is  computed. 

•  The  ratio  of  Q  to  is  computed,  and  a  random  variable  is  generated.  The 
pair  is  accepted  for  collision  \f  the  random  variable  is  less  than  that  ratio,  (If 
the  pair  is  not  accepted,  then  another  random  pair  is  selected.  The  process 
continues  until  a  pair  is  accepted.) 

•  For  an  accepted  pair,  the  collision  mechanics  are  computed  as  described  in 
Section  8. 

•  The  global  time  counter  is  incremented  by  VV^A/ ,  where  dJg  is  given  in 
Equation  (66),  and  is  the  lesser  of  the  two  weighting  factors. 

•  The  process  continues  until  the  global  time  counter  goes  beyond  the  overall  flow 
time.  At  that  point  the  collision  sampling  is  commenced  in  the  next  cell. 

•  When  all  celts  have  had  collisions  simulated,  then  the  code  proceeds  to  the 
translation  portion.  (See  Sections  I  and  6.) 
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7.5  Deviations  from  the  General  Procedure 


There  are  some  exceptions  to  the  above  relations  which  have  been  added  to  SOCRATES  in 
order  make  it  more  efficient.  These  exceptions  are  described  in  the  following  subsections. 


7.5.1  Cell  Specific  At„ 

Before  collisions  are  simulated  in  a  cell,  the  mean  residence  time  of  molecules  in  the  cell  is 
estimated  using  the  cell  dimensions  and  the  molecular  velocities.  When  collisions  are  simulated  in 
the  cell,  it  is  done  for  an  increment  of  the  global  time  counter  that  is  20%  of  this  mean  residence 
time  (but  no  less  than  At,J.  Collisions  are  not  again  simulated  in  the  cell  until  the  overall  flow 
time  has  caught  up  to  the  global  time  counter  for  the  cell. 

The  major  reason  for  doing  this  is  to  recognize  that  some  cells  will  tend  to  have  their 
molecules  remain  in  them  much  longer  than  others.  Cells  which  have  longer  molecular  residence 
times  will  tend  to  have  molecules  which  experience  more  collisions  within  the  cell.  When  the 
number  of  collisions  per  molecule  becomes  sufficiently  large,  it  can  be  assumed  that  the  molecules 
in  the  ceil  equilibrate  with  each  other,  and  the  equilibrium  sampling  procedures  described  in 
Section  9  can  be  applied.  Since  these  relations  are  much  faster  than  direct  collision  sampling,  it  is 
hishiy  desirable  to  apply  them  whenever  they  are  valid. 

For  unsteady  simulations  the  ceil  specific  At,„  is  not  applied  since  it  might  result  in  a 
temporal  blurring  of  the  solution. 


7.5.2  Relaxation  of  Q...'. 

The  current  value  of  Q,,^,  in  a  cell  Is  reduced  by  a  factor  of  0.95  if  20  or  more  potential 
collision  pairs  are  rejected  in  a  row.  ITie  rejection.'  of  coiitsion  pairs  can  become  the  most  time 
consuming  part  of  the  simulation,  and  a  large  value  of  *3xacefbates  the  problem.  This  change 
means  that  a  celt  is  nut  permanently  penalized  for  a  single  event  that  once  occurred  in  it,  but  the 
change  in  is  not  so  great  u.s  to  invalidate  the  pair  selection  probability.  Iliis  moditlcatiun 
can,  under  some  collision  dominated  circumstances,  result  in  an  order  of  magnitude  increase  in 
computational  speed. 


7.5.3  Maximum  Time  Counter  Increment 


Since  Atg  Is  inversely  proportional  to  relative  velocity  (Equation  (66)),  when  a  very  low 
velocity  collision  dues  occur,  it  can  result  in  very  large  increment  to  the  collision  time  counter, 
which  effectively  turns  off  collisions  in  the  cell  for  a  long  time.  Although  this  is  statistically 
proper  in  the  long  run,  it  can  result  in  a  substantial  statistical  fluctuation  in  the  short  run.  llie 
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codes  do  not  allow  a  collision  time  increment  to  be  greater  that  At,„,  the  overall  step  that  is  used  in 
the  solution. 


The  limitation  on  Atg  is  achieved  by  decreasing  the  weighting  factor  of  the  collision  below 
the  weighting  factor  of  either  of  the  two  colliding  molecules.  The  maximum  collision  weighting 
factor,  (WJ^,  is  given  by 


=  yAt^Vn^ffijC, 


(68) 


The  weighting  factor  that  is  applied  to  a  collision  is  actually,  therefore,  the  smaller  of  {Wj,  Wj, 
(^c)m«x)-  represents  this  value,  then  a  collision  counts  as  W^.  events.  In  order  to  maintain 
an  overall  correct  simulation,  the  Q  described  in  Equation  (67)  is  actually  given  by  the  relation 


^  “  W,. 


(69) 


and  the  lime  increment  applied  to  the  global  time  counter  is  W^.Atg.  (1'he  two  molecules  then  have 
their  state  vectors  updated  as  a  result  of  the  collision  with  probabilities  of  Wg/W;  and  W^/Wj 
respectively.)  .Most  of  the  time  is  equal  to  W^,  and  this  procedure  reduces  to  that  given  above; 
however,  the  problem  of  occasional  large  time  increments  Is  eliminated. 


7.5.4  Separation  of  Major  and  Minor  Species 

A  problem  arises  when  a  cell  happens  to  contain  a  single  molecule  of  one  mqjor  species  and 
all  other  molecules  in  the  cell  are  minor  species  with  weighting  factors  considerably  less  than  that 
of  the  first  molecule.  (In  treating  minor  species,  the  ratio  of  weighting  factors  may  be  as  much  as 
1000  or  more.)  Since  a  molecule  cannot  collide  with  itself,  collisions  between  major  species  can 
not  occur  in  such  a  cell.  The  result  is  that  the  contribution  of  major  species  collisions  to  the 
overall  time  counter  are  unobtainable;  and  the  entire  collision  time  increment  has  to  be  made  up 
with  collisions  between  the  single  major  species  and  one  or  other  of  the  minor  species.  (Collisions 
between  minor  species  were  rare  since  pairs  are  selected  with  a  probability  which  is  proportional  to 
the  greater  weighting  factor  -  see  above.)  The  result  of  this  problem  is  that  vastly  too  many 
collisions  are  simulated  between  the  major  and  minor  species.  This  is  both  unphysical  and 
numerically  inefficient. 

The  solution  to  the  problem  is  twofold;  I)  Logic  is  in  the  culiision  routines  to  recognize 
when  this  problem  occurs;  and  2)  ITie  global  time  counter  is  rcdetlned  in  such  situations.  TItc 
different  time  counter  is  achieved  simply  choosing  a  differem  dcfmiiion  for  the  D,j  coefficients 
appearing  in  Equation  (61).  Rather  than  taking  a  weighted  average  over  uii  collision  classes,  it  is 
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possible  to  take  a  weighted  average  over  just  those  collision  classes  which  involve  a  collision 
between  a  major  and  a  minor  species.  If  n|  represents  the  total  major  species  number  density  and 
nj  represents  the  total  minor  species  number  density,  then  Dy  is  defined  for  this  case  to  be 

Dy  =  njUj  ,  (70) 

rather  than  the  value  given  in  Equation  (63).  The  implied  increment  fur  the  global  time  counter  is 
given  by 


At.  =  -r. — -  . 

*  Vn,nja|jC, 


(71) 


This  counter  is  then  only  applied  for  collisions  between  the  species  declared  to  be  major  and  the 
species  declared  to  be  minor,  but  the  Increment  that  is  applied  is  much  larger  than  if  the  weighting 
were  over  ail  collision  classes.  Note  that  collisions  between  minor  species  still  occur  -  they  just 
don't  affect  the  collision  time  counter,  it  should  be  stressed  that  this  modification  is  only  applied 
for  the  special  case  of  a  cell  that  has  a  single  major  species  molecule  (defined  as  a  weighting  factor 
at  least  ten  times  greater  than  that  of  the  other  molecules).  The  result  of  this  modification  is  the 
return  of  the  correct  collision  frequency  to  the  simulation  and  the  removal  of  a  substantial 
numerical  problem. 
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8.  COLLISION  MECHANICS 


8. 1  Relations  for  Elastic  Collisions 

The  purpose  of  this  section  is  to  present  relations  appropriate  to  the  simulation  of  a  collision 
in  the  SOCRATES  code.  (The  question  of  how  molecules  are  selected  for  collisions,  which  is 
crucial  to  the  proper  simulation  of  collision  frequency,  was  taken  up  in  the  previous  section.) 
Conservation  of  momentum  implies  that  the  center-of-mass  velocity  of  the  collision  pair  is 
unchanged  by  the  collision;  and  conservation  of  energy  then  implies  that  the  magnitude  of  the 
relative  velocity  between  the  collision  partners  is  also  unchanged  by  the  collision.®  Since  the 
collision  is  treated  as  a  statistical  event,  all  that  remains  is  to  select  the  direction  of  the 
post-collision  relative  velocity  vector  from  the  correct  distribution.  As  mentioned  in  Section  it, 
collisions  in  the  VHS  model  are  treated  as  hard  sphere  collisions  when  they  occur  (though  they  do 
not  occur  with  the  same  velocity  dependence  as  do  hard  sphere  collisions).  Hence,  as  far  as  the 
collision  mechanics  is  concerned,  the  model  is  a  hard  sphere  model.  For  hard  sphere  molecules,  all 
directions  for  the  post-collision  relative  velocity  vector  are  equally  likely.  This  is  the  chief 
computational  simplicity  of  the  VHS  model. 

Let  the  two  molecules  be  identified  by  subscripts  1  and  2,  with  m  and  7  denoting  their 
masses  and  velocities.  If  i  and  f  indicate  initial  and  final  states,  then  the  relations  for  the  collision 
can  be  summarized  via: 


**  mi  +  im  ’ 

(72) 

Vr  =  IVii-Vjil  , 

(73) 

cos(fi)  =  1  -  2|S  , 

(74) 

sin(fi)  =  Vl-cos2(e)  , 

(75) 

<l>  =  2t0  , 

(76) 

-  v,lcos(fl),  sin(^)cos((i>),  slti(5)sln(<^)J  , 

(77) 
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and 


-  .  "!2Vrf 

Vlf  ''cm  +  +  m2  ’ 


(78) 


-  _  -  miVrf 

^2{  ~  '^cm '  ni|  +  m2 


(79) 


Where,  again,  indicates  a  random  variable  which  is  evenly  distributed  on  the  interval  zero 
to  one.  Each  time  that  appears  a  distinct  evaluation  of  the  random  variable  is  implied. 


8.2  Effect  of  Coordinate  System 

Note  that  the  expression  for  the  post-collision  relative  velocity  vector  (Equation  (77))  is  not 
coordinate  system  specific.  The  indicated  vector  components  can  apply  to  any  locally  orthogonal 
coordinate  system,  since  the  direction  implied  is  random.  The  convenient  coordinate  system  to  use, 
of  course,  is  the  coordinate  system  used  to  define  the  velocity  elements  of  the  state  vector. 

Although  Equation  (77)  is  independent  of  coordinate  system,  there  is  a  source  of  error 
which  is  dependent  on  coordinate  system.  This  error  arises  from  a  basic  premise  of  the  direct 
simulation  Monte  Cario  method,  namely  that  position  in  the  ceil  is  ignored  when  seiecting  collision 
partners.  If  the  velocities  are  expressed  in  a  coordinate  system  which  has  spatially  varying  basis 
vectors,  then  differences  in  position  bttween  the  two  molecules  can  imply  an  erroneous  difference 
in  velocity. 

SCX^RATES  makes  use  of  the  effect  of  spatially  varying  basis  vectors  to  solve  an  otherwise 
difficult  problem.  The  problem  arises  due  to  the  presence  of  concentrated  sources  of  contaminants, 
such  as  thrusters  and  evaporator  vents,  which  are  modeled  as  point  sources  producing  molecules 
traveling  (initially)  directly  away  from  the  source.  Since  there  is  no  length  scale  to  a  point  source, 
the  assumption  that  properties  are  constant  for  the  cells  in  the  immediate  vicinity  of  the  source  must 
be  invalid.  This  can  result  in  improper  collision  sampling  if  special  care  is  nut  taken. 

If  the  velocities  are  expressed  in  Cartesian  coordinates,  for  instance,  then  two  molecules 
selected  ftom  different  positions  within  the  cell  containing  such  a  point  source  can  have  a 
substantial  relative  velocity.  This  relative  velocity  is  illusory,  however,  since  it  merely  results 
from  the  assumption  of  a  point  source  and  the  neglect  of  spatial  differences:  there  should  not  be 
collisions  based  on  this  relative  velocity  since  Uie  molecules  are,  in  fact,  heading  away  from  each 
other. 


A  simple  resolution  to  this  problem  is  to  express  the  source  velocity  elements  in  spherical 
polar  coordinates.  In  these  coordinates,  every  molecule  leaves  Ute  poi,nt  source  with  the  same 
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velocity  in  the  direction  of  the  spherical  radius  vector.  Expressed  in  spherical  polar  coordinates, 
the  relative  velocity  disappears.  SOCRATES  transforms  velocity  vectors  to  spherical  polar 
coordinates  for  cells  in  the  vicinity  of  point  sources  (specifically,  when  the  total  number  density  is 
greater  than  three  times  the  ambient  number  density).  Collisions  are  sampled  in  the  transformed 
coordinates,  and  then  the  velocity  elements  are  transformed  back  to  the  normal  representation  after 
collisions  have  been  sampled  for  the  cell  in  question. 


8.3  Simulation  of  Inelastic  Collisions 


SOCRATES  uses  the  Borgnakke  and  Larsen'*  phenomenological  model  for  transfer  of  energy 
between  internal  and  translational  modes.  In  this  model,  a  collision  is  assumed  to  be  either 
perfectly  elastic  or  perfectly  inelastic,  via  a  user  specified  probability.  A  perfectly  inela.stic 
collision  is  achieved  by  summing  the  total  pre-collision  energy  (internal  energy  of  both  molecules 
plus  the  translational  energy  of  their  relative  motion.  Equation  (5)),  and  then  assigning 
post-collision  values  from  the  equilibrium  distribution  for  collisions  with  that  total  amount  of 
energy,  taking  into  account  the  number  of  internal  degrees  of  freedom  in  the  two  molecules.  Note 
that  this  model  has  the  ability  to  relax  from  a  nonequilibrium  to  an  equilibrium  state  via  an 
effective  collision  number.  The  ability  to  exchange  internal  energy  in  such  a  manner  comprises  a 
significant  increase  in  capability  for  Monte  Carlo  codes  beyond  the  previous  models  where 
molecule.*  had  no  internal  energy.  It  is  this  capability  which  enables  the  codes  to  realistically 
predict  the  macroscopic  effects  of  polyatomic  gas  flow. 


Let  and  {2  the  number  of  internal  degrees  of  freedom  of  the  two  molecules  in  an 
inelastic  collision,  and  E,  be  the  total  collision  energy  defined  by 

E,  «  Eei  +  E,i  +  Eji  ,  (80) 

where  Ebi  is  the  initial  translational  collision  energy  defined  by  Equation  (S),  and  En  and  £2;  are 
the  pre-collision  internal  energies  of  the  two  molecules.  Using  the  procedures  present^  in 
Reference  5,  the  somewhat  cumbersome  expressions  given  in  Reference  4  can  be  recast  in  terms  of 
the  chi-square  distribution.  Post-collision  values  for  the  respective  energies  are  given  by 


X|E. 

X,  ■fX2-»-X3  ’ 


X2E, 

X,  +  X2-fX3  ' 


and 


(81) 

(82) 
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(83) 


Ecf  = 


X3E, 

X,  +  Xj  +  X3  ’ 


where  Xj  is  selected  from  a  chi-square  distribution  with  degrees  of  freedom,  Xo  is  selected  from 
a  chi-square  distribution  with  degrees  of  freedom,  and  X3  is  selected  from  a  chi-square 
distribution  with  2(2  - «)  degrees  of  freedom.  (Efficient  procedures  for  sampling  from  a  chi-square 
distribution  are  also  given  in  Reference  5.)  The  post-collision  translational  energy  is  then  used  to 
determine  a  new  relative  velocity  between  the  two  molecules.  With  this  new  relative  velocity,  the 
previous  relations  for  determining  the  post-collision  velocity  elements  of  the  molecules  apply  for 
inelastic  collisions  as  well  as  for  elastic  collisions. 

The  fact  that  the  translational  energy  is  selected  from  a  distribution  with  2(2  - «)  rather  than 
the  3  degrees  of  freedom  that  might  be  expected  for  translational  energy  merits  some  explanation. 
It  is  due  to  the  fact  that  these  molecules  are  not  random  samples  from  the  gas,  but  rather  special 
molecules  owing  to  their  being  the  product  of  a  collision.  This  point  can  perhaps  best  be  seen  by 
considering  microscopic  reversibility,  where  the  inverse  collision  occurs  with  the  same  rate  in 
equilibrium.  For  this  reverse  process,  molecules  participating  in  it  are  not  all  equally  probable, 
since  those  with  greater  relative  velocities  are  more  likely  to  collide.  Hence,  the  number  of 
degrees  of  freedom  does  take  on  the  value  three  for  the  special  case  of  u  equal  to  1/2,  which  is 
precisely  the  case  of  collision  frequency  being  independent  of  relative  velocity.  Translational 
energy  in  collisions  behaves  as  if  it  has  2(2  -  w)  degrees  of  freedom. 


8.4  Collisions  Between  Molecules  with  Distinct  Weighting  Factors 

There  is  an  obvious  problem  when  considering  a  collision  between  two  simulated  molecules 
with  distinct  weighting  factors,  since  they  represent  a  different  number  of  real  molecules.  If  Wy 
and  Wy  represent  the  weighting  factors  for  the  two  molecules,  with  Wy  being  greater  than  Wy, 
then  the  collision  is  generally  counted  as  Wy  "events".  (More  precisely,  the  weighting  factor 
applied  to  the  collision  is  generally  taken  to  be  Wy.)  This  is  accomplished  by  always  assigning 
post-collision  velocity  and  energy  components  to  the  state  vector  of  the  molecule  with  the  smaller 
weighting  factor,  but  only  changing  the  components  of  the  molecule  with  the  greater  weighting 
factor  some  of  the  time.  The  probability  that  the  molecule  with  the  greater  weighting  factor  will 
have  its  components  changed  is  simply  Wy/Wy.  Statistically,  this  means  that  fur  a  large  number  of 
simulated  collisions,  each  such  simulated  collision  will  average  out  to  Wy  real  collisions  for  each 
species,  even  though  their  weighting  factors  differ.  It  should  be  noted  that  this  does  violate 
conservation  of  momentum  and  energy  on  an  individual  collision  basis,  but  those  quantities  are 
conserved  in  the  aggregate  over  a  large  number  of  collisions. 

In  some  cases  the  collision  is  assigned  a  weighting  factor  W^  which  is  less  than  either  of  Wy 
or  Wy.  When  this  is  done,  the  velocity  components  and  internal  energies  of  the  two  molecules  are 
changed  with  a  probability  of  W^/Wy  and  W^/Wy,  respectively.  (See  Section  7.) 
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8.5  Reactive  Collisions 


Reactive  collisions  can  be  simulated  directly.  The  treatment  of  reactive  collisions  is  similar 
to  that  for  inelastic  collisions,  except  that  a  heat  of  reaction  is  added  to  the  total  energy  expressed 
in  Equation  (80).  Reactive  collisions  can  result  in  the  disappearance  of  reactant  molecules,  with  the 
post-collision  state  being  applied  to  the  product  molecules. 


8.5.1  Types  of  Reactive  Collisions 

SOCRATES  has  a  fairly  comprehensive  chemistry  package  which  is  capable  of  handling  a 
variety  of  reactive  collisions.  Generally  speaking,  a  reactive  collision  is  an  event  which  occurs  due 
to  collisions  with  a  probability  that  depends  on  the  velocity  (or  energy)  of  the  collisions.  The 
following  generic  types  of  reactions  are  treatable: 

1.  Specific  Bimolecular  Reactions,  i.e.,  reactions  of  the  form 

A  -f  B  -*  C  +  D  , 

where  A,  B,  C,  and  D  are  particular  species.  An  example  of  a  reaction  of 
this  type  is 

O  +  HaO-^O  +  HjO*  . 

(In  this  example,  the  vibratlonally  excited  state  of  water,  HoO*,  is  treated  as 
a  distinct  species.) 

2.  Generic  Bimolecular  Reactions,  i.e.,  reactions  of  the  form 

A  +  M  -*  B  +  M  , 

where  A  and  B  are  particular  species,  and  M  can  be  any  species.  An 
example  of  a  reaction  of  this  type  is 

HoO  +  M-HjO*  +  M  , 

which  is  similar  to  the  previous  reaction  except  that  now  any  molecule  can 
serve  to  excite  the  water  molecule. 

3.  Dissociation  Reactions,  i.e.,  reactions  of  the  form 

A  +  M-»C  +  D  +  M  , 

where  M  is  any  molecule,  and  C  and  D  are  the  fragments  of  A  that  result 
from  dissociation.  An  example  of  a  reaction  of  this  type  is 

02  +  M-*0  +  04'M 
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8.5.2  Reactive  Collision  Probability 


The  Monte  Carlo  program  simulates  all  of  the  above  reaction  types  by  calculating  a  reactive 
cross  section  which  is  a  function  of  the  relative  collision  energy.  When  a  collision  occurs,  the 
reaction  is  simulated  with  a  probability  which  is  proportional  to  the  ratio  of  the  reactive  to  collision 
cross  section  at  the  relative  velocity  for  the  collision. 

There  are  two  options  for  specifying  the  reactive  cross  section.  The  first  is  to  specify  an 
Arrhenius  rate  constant,  of  the  form 

k,  =  A'Texpf-E./RoT)  ,  (84) 

where  E,  is  the  activation  energy  of  the  reaction  and  A  and  n  are  parameters  of  the  relation.  (R^  is 
the  gas  constant  and  T  is  temperature.)  The  unique  reactive  cross  section,  cr*,  corresponding  to 
Equation  (84)  is  given  by^ 


* 


VrU 


V^(l  +  5y)A 
2Ro''r(n  +  3/2) 


(Ec  ■  E«)" 


> 


(85) 


where  dy  is  unity  for  like  reactants  and  zero  for  unlike  reactants,  F  represents  the  gamma  function 
and  Eq  is  the  collision  energy  given  by  Equation  (5).  Note  that  a  rate  constant  is  defined  in  terms 
of  an  equilibrium  velocity  distribution,  so  the  correspondence  between  Equations  (84)  and  (85)  can 
be  made.  There  is  no  requirement,  of  course,  that  the  reactive  cross  section  given  by  Equation  (85) 
be  used  only  in  equilibrium  situations.  When  this  option  is  used,  only  the  arrhenius  parameters 
A,n  and  E^^  need  be  specified;  the  program  automatically  computes  the  corresponding  reactive  cross 
section. 

For  some  reactions,  the  form  of  Equation  (85)  is  too  restrictive,  and  it  is  then  possible  to 
input  a  table  giving  the  reaction  cross  section.  The  form  of  the  table  is  of  the  same  functional  form 
as  Equation  (85),  namely  the  product  of  the  relative  velocity  times  the  reactive  cross  section  is 
given  as  a  function  of  relative  collision  energy.  Although  this  form  is  nut  standard,  it  is  far  more 
convenient  for  reactions  where  one  of  the  reactants  is  generic  ("M"),  since  there  is  no 
correspondence  between  collision  velocity  and  collision  energy  until  the  masses  of  both  reactants 
are  specified. 


8.5.3  Options  for  Simulating  Reactive  Collisions 

SOCRATES  has  distinct  options  for  simulating  reactive  collisions  which  are  reflective  of 
diffierent  anticipated  user  needs.  In  all  options,  the  sampling  of  the  reaction  rate  (if  it  is  being 
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performed)  is  done  the  same  way.  Whenever  a  collision  occurs  between  the  two  reactants,  the 
reactive  cross  section  is  calculated,  and  the  reaction  is  counted  with  a  weighting  factor,  W^,  given 
by 


W,  =  W, 


V,(T 


(86) 


where  Wg  is  the  collision  weighting  factor  (see  the  previous  section).  Hence,  even  though  the 
reactive  cross  section  may  be  significantly  smaller  than  the  collision  cross  section,  the  statistics  on 
the  reaction  rate  are  similar.  (The  statistics  for  the  reaction  rate  may  converge  slower  due  to  the 
velocity  dependence  of  the  reaction  cross  section;  but  not  due  to  its  absolute  magnitude.)  If  two 
molecules  can  participate  in  multiple  reactions,  statistics  are  kept  for  each  reaction. 

If  products  are  introduced  as  a  result  of  the  reaction,  they  can  be  introduced  at  every 
simulated  reactive  collision  with  a  weighting  factor  of  Wp  or  introduced  with  a  weighting  factor  of 
Wg,  but  only  W/Wg  of  the  time.  The  difference  depends  on  the  importance  of  tracing  product 
species  in  the  simulation.  The  former  approach  will  result  in  more  computational  effort  being 
spent  on  the  product  species,  but  it  will  give  better  statistics  on  them.  In  either  case,  reactants  are 
removed  from  the  simulation  with  a  probability  of  W/Wg  in  any  reactive  collision. 


9.  PROCEDURES  FOR  COLLISION  DOMINATED  FLOW 


One  of  the  major  difficulties  in  the  classical  Monte  Carlo  technique  is  the  attainment  of 
equilibrium,  where  the  collision  frequency  can  become  prohibitively  large  for  a  direct  simulation. 
There  are  two  basic  approaches  for  this  problem,  and  both  are  utilized  in  SOCRATES.  The 
equilibrium  modeling  is  only  applied  if  no  products  are  being  introduced  into  the  simulation  as  a 
result  of  chemical  reactions.  When  such  products  are  introduced,  collisions  are  always  sampled  for 
the  ftill  time  increment.  The  two  equilibrium  approaches  are  described  below. 


9. 1  Collision  Cutoff  Approach 

This  is  the  usual  method  of  dealing  with  a  collision  dominated  flow  field.  In  this  method, 
collisions  are  sampled  in  a  given  cell  only  until  enough  have  been  sampled  to  guarantee 
equilibrium.  Since  further  collisions  only  result  in  the  maintaining  of  equilibrium,  they  need  not 
be  simulated.  It  is  necessary,  of  course,  to  estimate  the  actual  collision  frequency  and  keep  proper 
count  of  the  collisions  (in  particular  the  excitations)  which  are  not  directly  simulated  in  order  to 
obtain  the  correct  collision  frequency.  Once  a  cell  has  had  its  collisions  cut  off  in  this  fashion,  a 
flag  is  set  for  it.  On  subsequent  calls,  the  equilibrium  aftermath  approach  (described  in  the  next 
section)  is  applied  to  the  cell.  (The  equilibrium  aftermath  approach  also  calculates  collision 
frequencies  and  switches  back  to  regular  collision  sampling  if  it  becomes  too  low.) 


9.2  Equilibrium  Aftermath  Approach 

It  is  possible  to  avoid  sampling  collisions  altogether  If  it  is  known  that  the  cell  is  in 
equilibrium.  This  is  done  by  calculating  the  total  cell  energy  and  momentum,  and  then  selecting 
post-coilision  velocities  from  the  appropriate  equilibrium  distribution.  Although  the  principle  is 
simple,  the  application  is  complicated  by  the  fact  that  molecules  in  the  cell  do  not  all  have  the  same 
statistical  weight.  In  some  ways  (e.g.,  the  determination  of  mean  velocity)  the  statistical  weight 
acts  like  an  effective  multiplication  of  molecular  mass  -  the  greater  a  molecule’s  statistical  weight, 
the  greater  its  contribution  to  the  mean  flow  velocity.  In  other  ways  (e.g.,  the  assignment  of 
post-coliision  thermal  velocities)  the  statistical  weight  does  not  affect  the  result  •  a  light  molecule 
should  generally  have  a  large  thermal  velocity  irrespective  of  its  statistical  weight. 

The  second  difficulty  with  the  formulation  of  this  approach  is  the  necessity  for  constraining 
the  total  energy  and  momentum  to  match  the  pre-collision  values.  Hence,  it  is  not  proper  to 
calculate  the  initial  energy  and  then  simply  sample  from  Maxwellian  distributions  with  the  same 
mean  energy,  since  such  a  distribution  has  a  finite  probability  of  producing  a  molecule  with  any 
energy  -  and  the  net  result  would  be  an  unacceptable  divergence  in  the  cell  energy  and  momentum 
from  the  initial  values.  Both  of  these  problems  are  avoided  in  the  steps  enumerated  below. 
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The  method  is  implemented  by  calculating  the  total  momentum  and  energy  in  the  cell  and 
then  "peeling  off'  one  molecule  at  a  time  from  the  others.  The  internal  mode  energy  and  energy  of 
relative  motion  for  that  molecule  (relative  motion  with  respect  to  the  remaining  molecules)  are 
selected  from  equilibrium  distributions,  except  that  a  scale  factor  which  is  proportional  to 
temperature  is  temporarily  left  undetermined.  The  process  is  repeated  sequentially  until  all  energy 
modes  have  been  assigned  values,  and  then  the  overall  multiplicative  constant  is  chosen  to  match 
the  known  total  energy  of  the  system  (thus  determining  the  temperature). 

Each  molecule  is  then  assigned  velocity  components  which  are  consistent  with  the  known 
relative  velocity  between  it  and  the  remaining  molecules;  and  then  conservation  of  momentum 
determines  the  mean  velocity  of  the  remaining  molecules.  Again,  the  process  is  repeated 
sequentially  until  all  velocities  and  internal  energies  have  been  assigned. 


9.2.1  Conserved  Quantities 

The  total  energy  and  center-of-mass  velocity  are  directly  computed  via  the  following 
procedure: 

1)  The  following  sums  are  evaluated,  summing  over  all  the  simulated  molecules  in  the 


cell: 

Si  »  ,  (87) 

Sj  «  J^WimiUi  ,  (88) 

S3  »  ^WimiVi  ,  (89) 

$4  =  J^WimiWi  .  (90) 

S,  =  2]wimi(u?  +  V?  +  w?)  ,  (91) 

and 

S,  -  ,  (92) 
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where  Wj,  nij  and  Ejj  are  the  statistical  weighting  factor,  the  mass  and  the  internal 
energy,  respectively,  of  the  ith  molecule;  and  Uj,  V;,  and  W;  are  its  velocity- 
components. 


2)  The  center  of  mass  velocity  components:  u*,  v*,  and  w*,  are  computed  via: 

u*  =  Sj/S,  ,  (93) 

V*  =  S3/S,  (94) 

and 

w*  =  S4/S,  .  (95) 

3)  The  total  translational  energy  of  relative  motion  between  the  molecules,  can  be 
represented  by: 

=  {J^WimiKUi  -  uV  +  (v; .  vV  +  (Wj  -  w*)2l  ,  (96) 


although  it  is  mure  easily  evaluated  by  the  mathematically  equivalent  expression 


Etn. 


S2^  +  S3^  +  S4^ 

Si  * 


(97) 


4)  The  total  cell  energy  is  therefore  given  by 

Eiot  ®  Sft  +  Ejn,  . 


(98) 


9.2.2  Center-of-Mass  Velocity  Distribution 

Given  that  a  group  of  N  molecules  is  in  equilibrium,  it  is  possible  to  determine  the  form  of 
the  distribution  function  fur  their  mass  averaged  velocity,  taking  into  account  their  different 
statistical  weighting  factors.  This  relation  Is  most  easily  demonstrated  by  relating  the  Maxwellian 
velocity  distribution  to  the  normal  distribution  of  statistics  and  then  utilizing  a  basic  statistical 
theorem. 

A  variable,  r,  is  di.stributed  according  to  a  normal  distribution  If  its  probahllity  density 
function,  f(r),  is  given  by 
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f(r)  =  exp(-^)/(ffV5;)  , 

2a- 

where  (P‘  is  the  variance  of  the  distribution.  (The  distribution  has  been  selected  with  zero  mean 
since  the  effect  of  non-zero  means  does  not  influence  the  velocity  differences  which  are  the  goal  of 
this  exercise.)  A  basic  result  of  statistics  is  that  if  rj  is  selected  from  a  normal  distribution  with 
variance  a|  and  t2  is  selected  from  a  normal  distribution  with  variance  oi,  then  the  variable  r3 
defined  by 


rj  =  ar,  +  jSrj  ,  (100) 

will  follow  a  normal  distribution  with  variance  03  where 

=  0^07  +  /S-02  .  (101) 

If  it  is  recognized  that  a  normal  distribution  is  the  same  as  the  Maxwellian  distribution  for  a  single 
velocity  component,  then  this  result  implies  the  distribution  for  the  center-of-mass  velocity 
components  obtained  by  averaging  over  N  molecules  as  in  Equations  (93)  •  (95).  The  result  is  that 
this  mean  velocity  follows  a  Maxwellian  velocity  distribution  appropriate  to  a  "super"  molecule 
whose  mass,  m,,  is  given  by 

"»,=■(  2][Wimi)V  (  ]^Wfmi)  .  (102) 

t  _  I  !  ^  I 


(The  temperature  of  the  distribution,  of  course,  is  the  same  as  that  used  to  select  the  constituent 
molecular  velocities.)  Although  Equation  (102)  is  not  intuitively  obvious  (to  these  authors, 
anyway),  it  does  yield  some  expected  iimits.  If  all  of  the  weighting  factors  are  the  same,  then  m,  is 
the  sum  of  the  masses  of  the  individual  molecules.  However,  if  one  molecule’s  weighting  factor  is 
much  larger  than  the  others  (resulting  in  the  center-of-mass  velocity  of  the  group  being  essentially 
equal  to  that  molecule’s  velocity),  then  the  distribution  of  center-of-mass  velocity  is  the  same  as  the 
distribution  for  that  one  molecule. 


9,2.3  Molecular  Relative  Velocity  Distribution 

The  relative  velocity  between  an  individual  simulated  molecule  (referred  to  as  "molecule  J") 
of  mass  mj  witlt  respect  to  the  center-of-mass  velocity  of  N  other  molecules  wilt,  therefore,  have 
the  same  distribution  as  the  relative  velocity  between  that  molecule  and  another  molecule  of  mass 
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m,.  It  is  a  well  known  result  that  this  velocity  distribution  is  a  Maxwellian  distribution  appropriate 
to  a  molecule  with  a  reduced  mass,  /tj,,  given  by 


mm„ 
“  m+n^ 


(103) 


Put  in  terms  of  the  chi-square  distribution  which  is  used  extensively  in  the  Monte  Carlo  model,  u^, 
(the  square  of  the  relative  velocity  between  molecule  j  and  the  center-of-mass  velocity  of  the  other 
N  molecules),  can  be  expressed 


(104) 


where  is  a  variable  selected  from  a  chi-square  distribution  for  the  relevant  three  translational 
degrees  of  freedom,  [k  Is  Boltzmann’s  constant,  and  T  is  the  (as  yet  undetermined)  temperature.] 


9.2.4  Translational  Energy  of  Relative  Motion 

The  total  translational  energy  of  the  molecules  (which  can  be  expressed  as  in  Equations  (87) 
•  (96)  above,  summing  over  all  N+I  molecules)  can  be  algebraically  recast  in  a  form  which 
specifically  shows  the  contribution  of  relative  motion  between  molecule  j  and  the  N  other 
molecules.  Specifically, 


(105) 


In  Equation  (105),  is  the  translational  energy  which  would  result  if  only  the  N  molecules  were 
included  in  the  previous  sums,  and  E^.^|  is  the  value  obtained  with  all  of  the  molecules.  The 
factor  in  the  difference,  qj,  is  the  "reduced  weighted  mass"  between  molecule  j  and  all  of  the  other 
molecules,  i.e., 


’J  S,+W,mj  ■ 


(106) 


where  S|  is  as  defined  in  Equation  (87).  applying  the  sum  to  the  N  remaining  molecules,  it  is 
crucial  to  note  that  m,,  as  dellned  in  Equation  (102),  determine.N  the  distribution  of  relative 
velocities  between  mj  and  the  other  N  molecules;  but  ijj,  as  defined  above,  determines  the  amount 
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of  energy  associated  with  that  relative  velocity.  Combining  Equations  (104)  and  (105)  gives  the 
translational  energy  contribution,  E(j,  as 


E,J  =  <lkT)ix^  .  (107) 

Note  that  this  effectively  gives  a  weighting  factor  associated  with  the  translational  energy 
contribution  of  molecule)  of  ijj//4j,. 


9.2.5  Determination  of  Temperature 

The  internal  energy  associated  with  molecule  j,  E^,  can  be  represented  simply  by 


Ey  =  (|kT)WjXy  . 


(108) 


where  Xy  is  a  variable  selected  from  a  chi-square  distribution  with  the  number  of  degrees  of 
freedom  appropriate  to  molecule  j*s  internal  modes.  As  discussed  above,  this  process  Is  then 
repeated  sequentially  for  each  molecule  in  the  cell.  Note  that  "N*'  in  the  above  relations  refers  to 
all  remaining  molecules  which  have  not  had  their  energies  determined  yet.  This  means  that  m,,  for 
instance,  changes  with  each  molecule  since  It  Is  defined  via  a  sum  over  these  remaining  molecules. 
The  last  molecule  has  no  translational  energy  of  relative  motion  associated  with  it  since  there  are  no 
remaining  molecules  for  It  to  be  moving  with  respect  to.  ft  dues,  of  course,  have  internal  energy. 

Summing  all  of  the  Ey  and  Ey  and  equating  them  to  the  known  total  energy  E,o(  (as  given  In 
Equation  (98))  then  determines  the  temperature  of  the  system,  ft  is  noteworthy  that  this 
temperature  is  not  determined  Just  by  the  total  energy  of  the  sy.stem,  but  also  by  the  statistical 
sampling  process.  This  is  consistent  with  the  fact  that  any  temperature  cmtld  re.suli  in  the  particular 
observed  velocities;  although  some  temperatures  are  much  more  likely  to  produce  them  than  others. 

Once  the  temperature  is  defined,  then  the  sequential  relative  velocities  squan«l  u^, 
(Equation  (104))  and  internal  energies  (Equation  (108))  are  determined,  it  is  then  a  simple  matter 
to  go  bMk  and  apply  these  values  to  select  individual  molecular  velocity  components  via  the  same 
procedure  described  in  Section  8. 


9.3  The  Number  of  Collisions  Required  to  Achieve  Equilibrium 

The  number  of  collisions  requirt^  to  achieve  equilibrium  depends  on  the  model  being 
employed  and  the  criterion  for  equilibrium.  (Equilibrium  is  apprttached  asymptotically  and,  as 
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such,  could  be  regarded  as  an  ideal  limit  which  is  never  realized.)  The  model  being  employed,  as 
discussed  in  Subsections  3.4  and  8.3  is  that  of  Reference  4.  In  this  "statistical  collision"  model,  a 
fraction,  a,  of  the  collisions  are  taken  to  be  "perfectly  inelastic";  that  is,  in  such  collisions  all 
translational  and  rotational  energy  of  the  colliding  molecules  is  made  available  for  distribution  to 
the  post-collision  state  vectors,  taking  into  account  the  number  of  translational  and  internal  degrees 
of  freedom.  The  rest  of  the  collisions  are  taken  to  be  completely  elastic,  with  no  interchange 
taking  place  between  the  translational  and  rotational  energy  modes.  The  parameter  of  the  model, 
a,  should  be  chosen  to  match  available  data  for  rotational  relaxation. 

Within  the  context  of  this  model,  the  question  to  be  addressed  is  how  many  collisions  are 
required  in  a  cell  before  the  model  predicts  that  it  is  essentially  in  equilibrium.  The  question  can 
be  made  independent  of  a  if  it  is  phrased;  "How  many  inelastic  collisions  per  molecule  must  be 
simulated  before  the  cell  can  be  considered  to  be  in  equilibrium?".  This  question  is  suitable  fur 
direct  investigation  with  the  model,  and  a  test  calculatiun  was  performed  to  answer  it.  The  test 
calculation  indicated  that  the  equilibration  is  90%  complete  after  approximately  3.08  inelastic 
collisions  (on  the  average)  for  each  molecule.  This  seems  to  represent  a  reasonable  point  at  which 
to  say  that  further  collision  simulation  is  unnecessary,  although  the  cutoff  is  of  necessity  somewhat 
arbitrary.  This  number  of  inelastic  collisions  serves  as  a  useful  benchmark  in  the  comparison  of  the 
collision  cutoff  and  equilibrium  aftermath  approaches,  and  it  also  serves  to  define  when  the 
application  of  the  equilibrium  aftermath  tq)pruach  is  vaiid. 


9.4  Method  Comparison 

Test  runs  were  run  where  the  collision  cutoff  approach  was  utilised  for  3.08  inelastic 
collisions  per  molecule.  (Since  a  b  0.2  was  used,  this  corresponded  to  about  IS  total  collisions 
per  molecule,)  The  time  required  to  oimpute  the  relaxation  via  collisions  was  then  compared  to 
the  time  required  to  utilize  the  equilibrium  aftermath  approach.  The  result  was  that  the  equilibrium 
aftermath  approach  was  almost  an  order  of  magnitude  (a  factor  of  9)  faster  in  achieving  the  same 
result.  This  ratio  will  no  doubt  vary  with  computer  and  spectilc  calculation  being  performed;  but  it 
Is  highly  likely  that  the  equilibrium  aftermath  will  always  come  out  considerably  faster.  It  is  for 
this  reason  that  the  mdtod  was  implemented  in  SOCRATES. 
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10.  STATISTICAL  SAMPLING  OF  OUTPUT 


10. 1  General  Considerations 

It  is  safe  to  say  that  the  molecular  state  vectors  as  they  exist  in  the  computer  do  not 
comprise  the  usual  desired  output  of  the  procedure.  With  rare  exceptions,  it  is  the  macroscopic 
quantities  such  as  temperature,  density,  mean  flow  velocity,  etc.  which  are  of  interest  -  not  the 
microscopic  quantities  represented  by  the  state  vector  of  an  individual  simulated  molecule.  The 
generation  of  the  desired  output  requires  that  the  macroscopic  quantities  of  interest  be  represented 
in  terms  of  statistical  sums  of  the  available  microscopic  quantities;  and  it  is  the  main  purpose  of  this 
section  to  present  these  correspondences.  All  sums  are  kept  in  terms  of  "real"  molecules  and 
events,  i.e.,  the  current  weighting  factors  are  included  in  the  sums.  This  is  essential  since  the 
weighting  factor  determines  the  statistical  importance  of  a  given  molecule.  Since  the  weighting 
factors  are  dynamically  and  unpredictably  adjusted  as  the  solution  progresses,  it  would  not  be 
possible  to  go  back  and  add  in  the  effect  of  weighting  factors  a  posteriori. 

In  general,  it  must  be  decided  ahead  of  time  exactly  what  output  is  desired  from  the  code, 
and,  therefore,  what  statistical  sums  should  be  kept  to  generate  it.  There  is  a  vast  amount  of 
potential  information  in  the  simulation,  and  it  is  not  reasonable  to  store  ail  possibly  interesting 
quantities  in  all  runs.  On  the  other  hand,  it  is  wasteful  to  completely  rerun  a  case  Just  because  the 
user  decides  there  was  an  additional  quantity  he  was  interested  in.  The  selection  of  output  for  a 
given  run,  therefore,  unavoidably  requires  user  judgment.  Once  the  user  has  decided  upon  the 
required  output,  the  determination  of  which  statistical  sums  are  required  is  done  automatically  by 
the  code.  Care  is  taken  to  make  sure  that  a  statistical  sum  is  nut  duplicated  internaliy  if  it  is 
required  by  more  than  one  requested  output  quantity. 

Some  initial  words  of  caution  are  required.  By  its  nature,  the  direct  simulation  Monte  Carlo 
method  works  with  far  fewer  molecules  than  nature  does  and  it,  therefore,  exhibits  considerably 
greater  statistical  variation  in  its  macroscopic  predictions.  To  reduce  these  variations,  the  code  is 
run  repeatedly  for  the  same  case,  increasing  the  statistical  base  from  which  the  macroscopic  output 
is  derived.  Useful  results  can  usually  be  obtained  with  a  modest  computational  effort.  However, 
this  statement  must  be  tempered  by  a  realization  of  the  convergence  rate  for  Monte  Carlo  sampling. 
Basically,  the  statistical  error  in  the  output  converges  as  one  over  the  square  root  of  the  sample  size 
(or  run  time).  Hence,  if  a  solution  looks  good,  but  the  user  decides  he  would  like  one  more 
significant  digit  (i.e.,  he  would  like  the  .statistical  error  to  be  reduced  to  0. 1  times  its  current  value) 
it  would  require  that  the  run  time  be  increa.sed  by  a  factor  of  1001  It  cun  be  seen  that  the  desire  for 
more  accuracy  can  quickly  turn  the  most  efficient  code  into  a  money  gobbling  nightmare.  When 
using  a  Monte  Carlo  technique,  one  must  accept  some  statistical  scatter  in  the  output. 
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10.2  Sampling  of  Instantaneous  Volumetric  Output  Quantities 

Instantaneous  volumetric  output  quantities,  such  as  density,  temperature  and  velocity,  can  be 
determined  by  examining  the  moiecuiar  state  vectors  at  a  particuiar  time  in  the  simuiation.  The 
code  pauses  in  the  simuiation  and  uses  the  moiecuiar  state  vector  eiements  to  add  vaiues  to  the 
statisticai  sums  appropriate  to  the  various  ceils  and  the  particular  time  that  it  paused.  It  then 
proceeds  with  the  simuiation  until  the  next  sampling  time.  As  the  code  goes  through  its  successive 
runs,  it  stops  at  the  same  points  in  the  simulation  every  time  and  adds  to  the  statistical  base  for  the 
sums.  (For  steady  state  cases,  it  simply  does  it  repeatedly  after  the  initial  transient  has  died  down.) 
The  items  listed  below,  with  their  statistical  definitions,  are  selectable  as  output  requests  in 
SOCRATES.  Summations  are  performed  over  all  applicable  simulated  molecules,  which  include 
separate  runs. 

•  TOTAL  NUMBER  DENSITY 

n  =  vl'-  <l(») 

•  MEAN  MOLECULAR  WEIGHT 

m  «  (IlO) 

5| 

•  X  VELOCITY  COMPONENT 


8 


h 

Se 


(IN) 


•  y  VELOCITY  COMPONENT 

S4 

Vy  W  (112) 

•  z  VELOCITY  COMPONENT 

S) 

v»  «  ^  (113) 
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•  OVERALL  TRANSLATIONAL  TEMPERATURE 


1  8,2  +  S42  +  S52 

T  I  > 


•  TRANSLATIONALTEMPERATUREIN  jTH  DIRECTION 


I  5.2 


•  INTERNAL  MODE  TEMPERATURE 


T,  = 


2S9 

Ro^iO 


where  the  indicated  sums.  S|(,  are  defined  by; 

Si  =  . 

i 

i 

53  «  y^WimiV^  , 

54  *  J^WimiVyi  , 

i 

55  ■  JJWjmjVji  , 

i 

i 

S,  »  * 

i 

S|  ■  JJwjmjVji  . 

i 


(114) 

(115) 

(116) 

(117) 

(118) 

(119) 

(120) 

(121) 

(122) 

(123) 

(124) 
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(125) 


S.  =  2 

i 


(126) 


With  the  exception  of  Equation  (1 10),  all  of  the  above  quantities  can  also  be  defined  and  calculated 
for  any  specified  species.  The  sums  are  the  same  except  that  only  molecules  of  that  species  are 
considered.  Before  printing  output  quantities,  they  are  always  transformed  to  standard  dimensions 
firom  the  internal  dimensionless  variables. 


10.3  Sampling  of  Time  Averaged  Output  Quantities 

Some  additional  quantities  of  interest  are  not  sampled  at  a  separate  sampling  time  as 
described  above,  but  rather  as  the  simulation  evolves.  Examples  of  such  quantities  are  collision 
rates,  reaction  rates,  mean  velocities  between  molecules,  etc.  For  the  most  part,  these  quantities 
depend  on  the  relative  state  of  more  than  one  type  of  molecule,  and  they  are  by  their  nature 
expressed  as  average  values  over  a  finite  time  interval.  The  formulas  for  calculating  these 
quantities  are  no  more  than  event  counters,  and  will  not  be  included  here.  The  following  quantities 
are  currently  available  as  output: 

•  Mean  Relative  Velocity  Between  any  Two  Species; 

•  R.M.S.  Deviation  of  Mean  Relative  Velocity  Between  any  Two  Species; 

•  Mean  Product  of  Cross  Section  Times  Relative  Velocity  Between  any  Two 
Species; 

•  Collision  Rate  Between  any  Two  Species; 

•  Reaction  Rate  for  any  Chemical  Reaction; 

•  Reaction  Rate  for  any  Photochemical  Reaction; 

•  Flux  Rate  for  any  Species  on  any  Surface  Element. 

The  sampling  for  all  but  the  last  of  these  quantities  occurs  in  the  collision  simulation  routines.  As 
pairs  are  considered  as  possible  collision  partners,  statistics  are  kept  to  generate  the  first  three 
quantities.  Statistics  on  collisions  and  reactions  are  kept  as  they  occur,  and  the  last  quantity  is 
determined  in  the  molecule  advancemem  routines. 
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11.  INNER  SOLUTION  REGION  PROCEDURES 


11.1  Motivation  For  The  Inner  Solution 

The  separation  into  an  inner  and  outer  solution  is  necessitated  by  the  physical  and 
computational  realities  of  the  interaction  between  contaminant  molecules  and  the  surrounding 
atmosphere.  The  problem  is  created  by  the  combination  of  the  following  two  points:  1)  In  order 
to  describe  the  scattering  of  contaminant  molecules,  the  solution  region  must  extend  several  mean 
free  paths  from  the  source  of  the  contamination,  which  typically  means  tens  of  kilometers;  and  2) 
At  the  same  time,  in  order  to  describe  the  detailed  interaction  between  the  flow  field  and  the 
specific  geometry  of  the  vehicle  in  question,  cells  must  be  smaller  than  characteristic  length  scales 
of  the  vehicle.  Combining  these  two  requirements  in  a  three  dimensional  flow  situation  means  that 
many  more  cells  and  simulated  molecules  are  required  than  is  computationally  feasible.  However, 
the  two  requirements  naturally  lend  themselves  to  two  distinct  solutions,  which  are  referred  to  as 
the  "inner"  and  "outer"  solutions.  The  outer  solution  determines  the  scattering  of  the  contaminants 
by  the  atmosphere  and  doesn’t  require  solution  ceils  which  are  small  compared  to  the  vehicle,  since 
the  characteristic  length  scale  for  this  scattering  is  much  larger  than  the  vehicle.  Hence,  the  outer 
or  scattering  solution  can  be  carried  out  first,  with  alt  contaminant  sources  assumed  to  be  merely  at 
the  origin.  The  outer  solution  then  determines  the  rate  at  which  both  scattered  and  atmospheric 
molecules  will  be  approaching  the  vicinity  of  the  vehicle,  and  this  (rather  than  the  undistut'ocd  free 
stream)  can  then  be  used  as  a  boundary  condition  for  the  inner  solution,  allowing  its  outer 
boundaries  to  extend  to  Just  beyond  the  vehicle. 


11.2  General  Considerations 

The  boundary  conditions  for  solutions  via  the  Direct  Simulation  Monte  Carlo  method 
Involve  introducing  molecules  at  the  boundaries  which  have  .statistically  correct  fluxes,  velocity 
distributions,  and  internal  energies.  Only  the  inwardly  directed  portion  of  the  velocity  distribution 
function  is  relevant;  the  outward  portion  follows  naturally  from  the  solution. 

For  a  general  nonequilibrium  flow,  it  can  be  a  formidable  task  to  determine  the  complete 
inward  velocity  distribution  at  the  boundaries  so  that  molecules  may  be  introduced  properly. 
Generally,  therefore,  the  boundaries  are  merely  placed  far  enough  away  from  the  interaction  region 
so  that  the  flow  may  be  assumed  to  be  in  equilibrium  at  the  boundaries.  Since  an  equilibrium  flow 
has  a  velocity  distribution  which  Is  characterized  by  the  well  known  Maxwell-Boltzmann  function, 
it  is  then  straightforward  to  determine  a  statistically  correct  procedure  for  molecule  introduction. 

Such  a  procedure  is  clearly  not  adequate  for  the  inner  solution  described  above.  The 
ambient  flow  is  likely  to  be  disturbed  near  the  vehicle  where  the  inner  solution  boundary  is  to  be 
placed,  ami  the  flux  of  scattered  contaminants  back  Into  the  solution  region  would  be  completely 
lost  if  the  normal  boundary  condition  were  used.  The  procedure  that  is  used,  therefore.  Is  to  gather 
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relevant  statistics  in  the  outer  solution,  and  devise  a  boundary  condition  fur  the  inner  solution 
which  is  consistent  with  these  statistics. 


11.3  Velocity  Variations 

11.3.1  Assumed  Form  of  the  Velocity  Distribution  Function 

The  procedure  that  is  being  adopted  is  to  use  a  generalized  form  of  the  equilibrium 
distribution  function.  The  generalizations  that  are  applied  are  as  follows:  1)  Each  species  is 
allowed  to  have  its  own  mean  flow  velocity.  (In  equilibrium,  all  species  would  have  the  same 
mean  velocity.)  2)  Each  species  is  allowed  to  have  its  own  translational  temperature,  and  this 
temperature  is  allowed  to  be  distinct  for  the  three  coordinate  directions.  (In  equilibrium,  all 
species  would  have  the  same  temperature,  and  it  would  not  vary  with  coordinate  direction.)  3) 
Each  species  is  allowed  its  own  internal  mode  temperature.  (In  equilibrium,  this  temperature 
would  be  the  same  for  all  species,  equal  to  the  common  translational  temperature.)  Quantitatively, 
it  is  assumed  that  for  the  ith  species,  the  velocity  distribution  in  the  jth  direction  follows  the 
functional  form: 


fijCuy) «  («ij/v7)expl-a?(u^.u^|)2j  (127) 

where 

ay  «  Vmi/(2RoTy)  .  (128) 


In  these  relations,  m|  is  the  molecular  weight  of  the  ith  species  in  atomic  mass  units,  Rq  is  the 
universal  gas  constant,  and  Ty  and  TTy  represent  the  translational  temperature  and  mean  velocity, 
respectively,  of  the  ith  species  in  the  Jth  direction.  By  definition  of  fy,  the  probability  of  a 
molecule  of  species  i  having  a  velocity  component  in  the  jth  direction  between  uy  and  Uy+Auy  is 
fy(uy)Auy.  Under  this  assumption,  therefore,  the  equilibrium  boundary  condition  can  be  applied 
for  each  species  (including  those  not  present  in  the  free  stream)  by  determining  the  parameters  ITy 
and  Ty,  as  well  as  the  effective  density  of  the  species,  n^. 


11.3.2  Available  Statistics  From  the  Outer  Solution 

The  outer  solution  will  have  the  location  of  the  outer  boundary  of  the  inner  solution 
specified.  This  boundary  will  be  subdivided  into  rectangular  elements  which  may  (but  need  not 
necessarily)  be  particular  cell  boundaries  for  the  outer  solution.  For  each  rectangular  element, 
whenever  a  molecule  crosses  that  boundary  in  the  outer  solution,  moving  into  the  inner  solution 
region,  statistics  can  be  kept.  In  addition  to  Just  counting  the  molecules  of  each  species,  it  is  also 
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possible  to  keep  velocity  and  spatial  moments  of  the  flux  of  each  species.  This  section  will 
demonstrate  how  velocity  moments  can  be  used  to  determine  the  parameters  TTij  and  Ty,  as  well  as 
the  effective  number  density  nj,  appropriate  to  a  particular  rectangular  element.  Subsection  1 1 .4 
will  deal  with  spatial  variations  along  the  rectangular  element. 

Let  Tj  be  the  inward  normal  at  the  surface  element  and  To  and  T3  represent  two  other  unit 
vectors  such  that  T],  T2  and  T3  form  an  orthonormal  triple  of  unit  vectors.  Depending  on  the 
orientation  of  the  rectangular  element,  each  of  these  unit  vectors  will  be  either  aligned  or 
anti-aligned  with  one  of  the  three  basic  solution  unit  vectors:  Ty  and  Tj,.  Keeping  statistics  in 

the  T2  and  i  3  directions  is  straightforward,  following  the  procedures  described  in  the  previous 
section.  The  more  difficult  task  is  to  determine  the  effective  equilibrium  gas  quantities  for  the 
direction  aligned  with  the  normal  of  the  surface  element,  Tj.  As  molecules  cross  the  boundary 
segment,  the  following  statistics  (among  others)  can  be  kept  for  each  species: 

a)  qj;  The  flux  of  species  1  In  molecules/(cm^s). 

b)  <  Uii  > ;  The  mean  value  of  the  ist  velocity  component. 

c)  <  u^i  > ;  The  mean  square  value  of  the  Ist  velocity  component. 

Note  that  the  sampling  is  applied  only  for  molecules  which  cross  the  boundary  with  a  positive 
velocity  component  in  the  T (  direction,  so  that  all  molecules  included  in  these  statistics  will  by 
definition  have  a  positive  value  of  un.  Since  only  molecules  with  a  positive  U||  component  are 
included  in  the  statistics,  and  since  the  molecules  with  large  un  components  will  tend  to  flux  across 
the  boundary  faster  and  be  counted  with  a  larger  weight,  U  is  not  the  case  that  <Ufj>  is  equal  to 

“u* 


11.3.3  Representation  of  Available  Statistics 

For  a  Maxwellian  velocity  distribution  (Equation  (127)),  it  is  possible  to  directly  compute 
die  expected  values  of  the  available  statistics.  The  results  are; 


m 

%»ni|ui 


ijfu(utiWui,  , 


<»»it>  “  ^  Juftfu^lWUil 


(129) 


(130) 
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and 


00 

<“u>  = '^|«ufii(Uii)dUii  . 


(131) 


Substituting  the  expression  for  fj,  given  in  Equation  (127),  and  performing  a  little  algebra,  these 
expressioris  can  be  written: 


nili(-w) 

(132) 

.  s._ 

”  aiili(-w) 

(133) 

and 

(134) 

where  the  speed  ratio,  w,  is  given  by 

W  -  0(„Ui, 

(135) 

and  the  integrals  1|(  are  defined  by 

00 

W  ■  . 

b 

(136) 

The  first  few  1|,  are  directly  evaluated  to  give: 

fjCb)  *  *  bv7  erfc(b)l  . 

(137) 

^(b)  *  i(.2b  exp(-b2)  +  (1  +  erfc(b)l 

(138) 
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and 


13(b)  =  |[2(1  +b2)  exp(-b2) .  b(2b2+3)V^  erfc(b)]  . 


(139) 


11.3.4  Inversion  to  Determine  Equilibrium  Parameters 

Combining  Equations  (133)  and  (134),  it  is  possible  to  get  a  relation  that  involves  only  the 
sampled  variables  and  the  speed  ratio,  viz: 


r 


<U?l>/(<Ui,»2 


ll(-w)l3(-w) 

(I2(-W)r- 


(140) 


Since  r,  the  ratio  of  the  mean  square  velocity  to  the  square  of  the  mean  velocity  for  molecules 
fluxing  across  the  rectangular  segment,  is  a  function  only  of  w,  it  is  a  useful  place  to  start  the 
inversion.  As  can  be  seen  from  Figure  4,  r  is  monotonically  decreasing  with  increasing  w,  going 
from  an  asymptote  of  1 .5  for  w  -» -«  to  an  asymptote  of  1 .0  for  w  -►  +  00 .  Since  the  function  is 
monotonic,  it  can  be  inverted  unambiguously  as  long  as  r  falls  within  the  interval  (1.0  <  r  <  l.S). 
Although  the  functional  form  of  r(w)  is  somewhat  complex,  making  an  analytic  inversion 
impossible,  a  Newton-Raphson  iteration  starting  at  the  center  of  the  curve  converges  rapidly.  Since 
this  calculation  is  only  done  once  for  each  segment  to  determine  the  boundary  condition 
parameters,  the  time  required  for  this  inversion  is  negligible. 

Once  w  is  determined,  then  Equation  (133)  can  be  used  to  determine  the  corresponding 
value  of  0(|i,  which  defines  Tn  via  Equation  (128).  w  and  a;i  determine  TTn  via  Equation  (135), 
and  Equation  (132)  then  determines  n^. 

11.3,5  Possible  Errors 

In  some  cases,  the  procedure  described  above  will  not  work  directly,  since  the  sampled 
value  of  r  may  exceed  1.5.  This  can  be  demonstrated  by  considering  the  simple  case  of  two 
molecules,  with  velocities  U|  and  Uj.  The  value  of  r  sampled  In  this  case  approaches  2.0  as  either 
molecule's  velocity  becomes  much  larger  than  the  other.  Hence,  it  is  distinctly  possible  that  values 
of  r  in  excess  of  1 .5  will  be  encountered,  either  through  statistical  error  or  the  failure  of  the 
assumption  of  an  equilibrium  type  velocity  distribution.  (It  is  mathematically  impossible  to  get  a 
value  of  r  less  than  1 .0,  though  numerical  error  might  conceivable  cause  it  in  extreme  cases.) 

To  handle  these  cases,  as  well  as  deal  with  some  numerical  problems  that  occur  fur  large 
negative  w,  die  inversion  described  in  the  previous  section  is  performed  only  if  the  resulting  value 
of  w  is  within  the  range  (‘5  <  w  <  45).  If  r  is  such  that  w  would  full  outside  of  that  range,  w  is 
set  equal  to  ±5  (depending  on  whether  r  is  too  large  or  too  small),  and  the  inversion  procedure 
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SPEED  RATIO,  w 


Fi^e  4.  A  Graphical  Representatiun  ut'  the  Variable  r,  The  Ratio  of  the  Mean 
Square  Velocity  to  the  Square  of  the  Mean  Velocity  tbr  Molecules  from  a 
Maxwetlicn  Velocity  Distribution  Fluxing  Across  a  Plane,  as  a  Function  of 
the  Speed  Ratio,  w. 
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described  above  is  then  otherwise  followed.  The  effect  of  this  is  that  the  derived  boundary 
condition  will  match  the  statistically  observed  flux  and  mean  velocity,  but  will  no  longer  match  the 
mean  square  velocity. 


11.4  Spatial  Variations 

The  rectangular  element  described  in  Subsection  11.3.2  will  typically  be  a  cell  boundary  in 
the  outer  solution.  When  the  boundary  condition  is  applied  in  the  inner  solution,  it  may  comprise 
the  boundary  of  several  cells,  since  the  inner  solution  will  by  definition  have  a  finer  cell  mesh  than 
the  outer  solution.  It  would  be  possible,  of  course,  to  merely  apply  the  same  boundary  condition  to 
each  of  the  corresponding  inner  solution  region  cells,  but  with  a  little  additional  work  it  is  possible 
to  include  some  realistic  spatial  variation  in  the  rectangular  region,  which  may  be  used  to  provide 
adjacent  inner  region  cells  with  slightly  different  boundary  conditions.  The  procedure  to  do  this  is 
described  in  this  section.  This  procedure  allows  for  a  spatial  variation  in  flux  (and  therefore 
effective  number  density)  along  the  rectangular  region,  but  the  velocity  distribution  is  assumed  to 
be  spatially  constant. 


11.4.1  Reduction  to  Scaled  Variables 

Let  the  rectangle  of  interest  be  represented  by  the  general  variables  X  and  such  that 
(X|  <  ^  <  :i(2)  and  (t|  <  ^  <  ^2).  Then  it  is  convenient  to  introduce  the  scaled  variables,  X 
and  Y,  such  that 


and 


X  >  -1  +  2 


X-X| 

X2-X| 


Y 


■  -1  +2 


111 

Yj-t, 


(141) 


(142) 


The  scaled  variables  have  the  advantage  that  they  each  traverse  the  range  of  >1  to  -f  1 . 


11.4.2  Spatial  Moments 

As  molecules  traverse  the  rectangular  region  in  question,  it  is  possible  to  calculate  the  X  and 
Y  values  at  which  a  given  molecule  crosses  it.  Sums  of  these  values  can  also  be  ketu,  so  that  a 
statistical  determination  can  be  made  of  the  average  X  value.  <  X  > ,  the  average  Y  value,  <  Y  > , 
and  the  average  product  of  the  two.  <  XY  > .  At  the  same  time,  of  course,  the  average  flux  across 
the  rectangle  wilt  also  be  determined. 
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11.4.3  Analytic  Representation  of  Sampled  Moments 

Let  j(X,Y)  be  flux  distribution  in  the  rectangular  region  (molecules  per  unit  area  per  unit 
time).  The  number  of  molecules  crossing  the  rectangle  per  unit  time,  N,  is  therefore  given  by: 


dXdY  . 


(143) 


The  average  sampled  values  given  in  the  previous  section  can  therefore  be  represented  via 

(144) 

(145) 

(146) 

11.4.4  Assumed  Form  for  the  Flux  Distribution 

U  is  assumed  that  the  flux  distribution  can  be  adequately  represented  by  a  linear  distribution, 
viz. 

j(X,Y)  «  li|(ao  +  a,X  +  ajY  +  ajXY)  .  (147) 

When  this  form  is  substituted  imo  the  integt'als  deflned  in  Equations  (143)  •  (146).  there  results; 
*0-7  .  (148) 


I  I 

<X>  f  Jj(X.Y)XdXdY  , 
*''.1 .1 

I  1 

<Y>  Jj(X,Y)YdXdY 

"  .1 .1 


and 


<XY>  » 


"•I  -i 


XYdXdY 
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(149) 


3<X> 

a,  =  — 5— 

3<Y> 

and 

9<XY> 


(150) 

(151) 


Hence,  under  the  assumption  of  a  linear  flux  distribution,  the  spatial  moments  lead  immediately  to 
the  corresponding  linear  coefflcients. 


11.4.5  Possible  Problems 

Since  the  flux  that  is  discussed  above  is  the  one-way  flux  (i.e.,  only  counting  molecules 
traveling  in  one  direction  across  the  rectangular  segment),  it  is  by  deflnition  a  quantity  which  can 
never  be  negative.  However,  again  due  to  either  statisticai  error  or  the  invalidity  of  the  assumed 
fiinctional  form,  it  is  conceivable  that  the  distribution  derived  above  could  become  negative  at  some 
subregion  of  the  rectangle.  If  this  does  occur,  corrections  must  be  made  to  the  coefflcients. 

In  order  to  search  for  a  negative  subregion,  it  sufflces  to  investigate  the  corner  values,  since 
the  maximum  and  minimum  must  occur  on  the  corners  for  a  linear  flt.  The  corner  vaiues  are  given 
directly  by 


Jii  *  “  Si  •  82  +  Mj)  ,  (152) 

ji2  •  K*^*'*'  “  *^^*0  ‘  Uj  +  02  -  aj)  ,  (153) 

ill  ®  j(-*‘ l.*l)  “  N(“o  +  *1  *  ^2  •  aj)  (154) 


and 


J22  •*  1)  *»  l^(ao  +  a|  +  02  +  a,)  .  (155) 

If  a  negative  value  is  found,  then  use  can  be  made  of  tlte  inverse  rdaijoiiV.  which  glW  .^a 
finmion  of  the  corner  values,  i.e., 


(156) 


®  ■*•2  jai  122) 

4N 


*i  “  *112  + 121  + 122)  » 

4N 


*2  “  ’*■  -^12  ■  j2l  i22) 


and 


«3  “  *  jl2  ‘hi  +  122) 


(157) 

(158) 
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If  one  of  the  corner  values  is  negative,  then  it  can  be  artificially  increased  to  zero,  while  one  third 
the  antount  added  to  the  negative  corner  is  subtracted  from  each  of  the  other  corner  values.  This 
procedure  keeps  ag,  and  therefore  the  integrated  lota)  tlux,  consiant.  Once  the  corner  values  have 
been  acyusted  so  that  none  are  negative,  Equations  (IS6)  •  (1S9)  give  the  corresponding  linear 
coefficients. 


*60- 


12.  COMPARISON  OF  SHUTTLE  ENGINE  FIRINGS  IN  SPACE  TO  CODE 
PREDICTIONS 


12. 1  Introduction 

This  section  describes  a  model/data  comparison  for  some  visible  data  taken  by  the  Air  Force 
Maui  Optical  Station  (AMOS)  of  870  lb.  shuttle  RCS  firings  at  altitudes  around  320  km.  The  data 
are  unique  in  that  they  were  obtained  for  three  angles  of  attack;  0°,  90'“,  and  180°.  The  sequence  of 
events  was  to  fire  in  one  direction  for  three  seconds  and  then  shut  down  for  five  seconds  before 
firing  another  engine  in  another  direction.  The  data  which  are  being  utilized  in  this  paper  involve 
total  visible  plume  intensity  as  measured  by  an  S-20  photocati.ode,  and  its  spectral  dependence  as 
measured  by  a  spectrograph.  In  particular,  we  are  concentrating  on  a  peak  in  the  spectrum  at  6300 
A  which  is  being  attributed  to  the  forbidden  transition  between  the  0(‘D)  excited  electronic  state 
and  the  0{T)  ground  state.  This  emission  is  intriguing  because  the  long  (*  194  seconds^) 
radiative  lifetime  of  the  excited  state  would  normally  be  associated  with  a  very  weak,  spatially 
diffuse  signature.  However,  the  quenching  of  the  0(‘D)  state  is  sufficiently  rapid  at  this  altitude 
that  the  net  lifetime  of  the  excited  state  is  more  in  the  vicinity  of  one  second,  providing  a  much  less 
diffuse  signature  than  would  exist  in  the  absence  of  quenching.  Furthermore,  the  abundance  of 
atomic  oxygen  is  sufficient  to  produce  i  measurable  signature  even  at  this  reduced  lifetime. 

The  existence  of  a  single  set  of  well  characterized  unsteady  data  as  a  function  of  angle  of 
attack  provides  an  excellent  opportunity  to  do  model  comparisons  with  the  hope  of  discerning  the 
underlying  mechanism(s).  Two  mechanisms  were  hypothesized  as  potentially  being  responsible  for 
this  emission:  the  charge  exchange  reaction  between  atmospheric  O'*’  ions  and  H^O  from  the 
exhaust,  and  the  direct  collisional  excitation  of  the  0{'D)  state.  These  data  provide  a  .stressing 
validation  case  for  the  SOCRATES  code,  as  well  as  an  opportunity  to  learn  about  a  new  visible 
plume  emission  source. 


12.2  Experiment 

The  experiment  consisted  of  firing  specific  870  lb.  PRCS  liquid-fuel  (monomethyl 
hydrazine-N204)  engines  while  the  shuttle  orbiter  passed  over  AMOS  atop  Mt.  Haleakala,  Maui, 
HI  (21  °N,  204°E,  3000  m  altitude).  The  engines  were  fired  for  angles  of  attack  of  180°,  90°,  and  0° 
with  respect  to  the  ambient,  which  will  be  referred  to  as  ram.  perpendicular,  and  wake  burns, 
respectively.  The  shuttle  was  in  total  darkness  (the  solar  and  lunar  depression  angles  at  the  shuttle 
were  30°  and  59°,  respectively),  Imagery  data  were  obtained  using  8"  and  22"  acquisition 
tele.scopes.  Visible  spectral  data  were  obtained  using  a  spectrograph  with  its  own  foreoptics 
mounted  on  the  main  telescope.  The  spectral,  spatial,  and  temporal  properties  of  the  emissions 
were  measured  as  a  function  of  angle  of  attack  of  the  exhaust  products  by  the  atmosphere.  The 
total  plume  intensity  fur  an  S-20  photocathode  was  estimated  using  the  6°  GEODSS  telescope  which 
had  a  fixed  field  of  view.  This  field  of  view  was  sufficient  to  encompass  the  entire  plume,  and  its 
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fixed  nature  enabled  known  stars  to  be  used  for  calibration.  The  experiment  is  described  in  more 
detail  in  Reference  10. 


12.3  Data  Reduction 

On  thf  basis  of  comparing  against  a  known  star  in  the  GEODSS  field  of  view,  the  peak  total 
plume  emission  for  an  S-20  photocathode  was  estimated  to  be  150,  80,  and  12  W/sr  for  the  ram, 
perpendicular,  and  wake  burns,  respectively.  An  average  spectrum  shows  a  substantial  NH  peak 
around  3360  A,  a  peak  around  6800  A  due  to  second  order  diffraction  of  the  NH(A-*X)  fundamental, 
as  well  as  several  other  peaks.  For  the  present  purposes,  the  quantity  of  interest  is  the  total 
emission  from  the  6300  A  atomic  oxygen  line.  To  estimate  this,  the  spectrum  was  multiplied  by  an 
S'20  response  curve  and  integrated  over  the  spectral  range.  The  ratio  of  the  integral  between  6100 
and  6S00  A  to  the  total  integral  was  estimated  to  be  the  fraction  of  the  total  S-20  response  which  was 
due  to  the  6300  A  transition.  For  the  present  case  this  ratio  came  out  to  be  12  % .  It  should  be  noted 
that  a  more  proper  approach  would  be  to  do  the  same  averaging  for  spectra  resulting  from  burns  in 
the  three  directions,  since  the  varying  available  kinetic  energy  may  shift  the  relative  magnitudes  of 
the  emission  peaks  as  a  function  of  angle  of  attack. 

The  atomic  oxygen  Green  line  {0('S)-*0('D))  at  5577  A  Is  nut  observed  in  any  of  the  burns. 
Based  on  the  increased  activation  energy  required  for  excitation  of  the  0('S)  state  this  seems 
reasonable.  11»e  0{'S)  lies  at  4.2  eV  above  the  0(*P)  ground  state  as  compared  to  2.0  eV  for  the 
0(>D)  state.  For  the  0  and  90  degree  burns  the  maximum  available  collisional  energy  is  not 
sufficient  to  reach  the  0('S)  level.  At  1 80  degrees  there  is  enough  collisional  energy  to  excite  the 
Green  line  emission,  but  it  is  stilt  e.xpected  to  be  less  than  the  Red  line  signal.  The  issue  of  the 
relative  Green  line  intensity  for  the  retro-fire  case  is  planned  for  a  future  investigation. 


12.4  Modeling  Approach 

The  atmosplieric  and  exhaust  mole  fractions  of  the  species  considered  in  the  calculations  are 
givan  in  Table  2.  The  modeling  efforts  focused  on  simulating  the  overall  character  of  the  plumes 
in  the  three  orientations,  and  in  trying  to  understand  the  6300  A  data  feature,  which  is  being 
Surihuted  to  the  forbidden  transition  between  the  0(‘D)  excited  state  and  the  0(*P)  ground  slate. 
Since  the  0('D)  state  has  a  long  radiative  lifetime,  it  was  expected  that  quenching  would  be  an 
essential  element  of  the  calculations.  The  following  rate  constants  (in  cm^/s)  from  the  literature 
were  used: 


0(‘D)  -b  Nj  -*  0  +  Nj 

k“  »  2.3x10’“ 

(160) 

0(‘D)+ 0-0  +  0 

k‘-  =  8.8xl0‘- 

(161) 

0('D)  +  H20-20H 

k'^  »  2.2x10’*° 

(162) 
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Table  2.  Concentrations  of  Species  Carried  in  the  Model  Calculations. 


Species 

Atmospheric 

Mole  Fraction 

Exhaust 

Mole  Fraction 

0 

0.9076 

0.0 

N2 

0.0919 

0.309 

^  H2 

0.0 

0.187 

H2O 

0.0 

0.332 

CO 

0.0 

0.136 

C02 

0.0 

0.036 

O'*- 

0.0005 

0.0 

Published  rate  constants  or  cross  sections  for  the  coliisional  excitation  of  the  0('D)  state  could  not 
be  found,  so  microscopic  reversibility  was  invoked  to  estimate  the  tbliowing  excitation  rate 
constants: 


0  +  N2-»0('D)+  Nj 

k  =  I.28xl0-“exp(-^^) 

(163) 

0  +  0-0('D)  +  0 

k  =  4.9xio*'2exp(.~|~) 

(164) 

Reaction  (163)  was  also  assumed  to  apply  for  the  other  major  plume  constituent.  H2O.  The  charge 
excha;)ge  mechanism  for  producing  0('0)  was  represented  by 

0^  +  HjO-OCD)  +  HjO*  k  *  3.23xlO  ’T®-^exp(-i^)  .  (165) 

The  rale  constant  fur  Reaction  (165)  was  taken  from  measurtnl  values'^  with  the  activation  energy 
adjusted  to  reflect  the  additional  energy  required  ti»  create  the  0('D)  rather  than  the  0(*P) 
product. 


12.5  Results  And  Conclusions 

Grayscale  plots  are  given  in  Figures  5-7  showing  the  calculated  evolution  of  the  6300  A  ram 
plume  radiance  as  a  function  of  time.  There  is  a  substantial  qualitative  resemblance  to  the  data  in 
these  plots,  though  a  quantitative  comparison  of  these  images  is  difficult.  Recall  that  the  engine 
shuts  off  at  3.0  seconds,  so  the  spreading  and  reducing  of  the  signature  at  4.0  seconds  in  Figure  7 
is  to  be  expected.  Thu  total  calculated  6300  A  radiant  intensity  as  a  function  of  time  is  given  in 
Figure  8  for  the  three  angles  of  attack.  If  the  peak  values  from  this  figure  are  compared  to  the 
deduced  data  values,  then  the  result  is  the  comparison  shown  in  Figure  9.  The  agreement  in  Figure 
9  is  probably  fortuitously  gtxxl,  since  it  is  doubtful  that  even  the  input  quantities  were  known  that 
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Figure  S.  A  Grayscale  Plot  Showing  the  Calculated  6300  A  Plume  Radiance  for  a  Ram 
Bum  at  0.4  Seconds.  The  Shuttle  is  Located  at  the  Center  of  the  Plot,  and 
the  Wind  is  Approaching  from  the  Right. 


180°  R.O.fl.  AT  1.6  SECONDS 
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Grayscale  Plot  Showing  the  Calculated  6300  A  Plume  Radiance  for  a  Ram 
im  at  1.6  Seconds.  The  Shunie  is  Located  at  the  Center  of  the  Plot,  and 
e  Wind  is  Approaching  from  the  Right. 
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Fipire  7.  A  Grayscale  Plot  Showing  the  Calculated  6300  A  Plume  Radiance  for  a  Ram 
Bum  at  4.0  Seconds.  The  Shuttle  is  Located  at  die  Center  of  the  Plot,  and 
the  Wind  is  Approaching  hrom  die  Right. 
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Figure  8.  The  Calcutated  Total  6300  A  Radiant  Intensity  as  a  Function  of  Time  for  the 
Three  Bum  Directions. 
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Figure  9.  A  Conq>arison  Between  the  Calculated  and  Deduced  Peak  Total  6300  A 
Plume  Emission  as  a  Function  of  Angle  of  Attack. 


accurately.  Nevertheless,  the  agreement  is  impressive  and  strongly  supportive  of  the  hypothesized 
mechanisms. 

Given  confidence  in  the  basic  description,  it  is  then  possible  to  use  the  model  to  get  more 
quantitative  and  qualitative  understanding  of  the  underlying  physics.  The  major  difference  to  be 
expected  in  the  excitation  mechanisms  is  the  activation  energy  required;  the  charge  exchange 
mechanism  requires  a  lower  activation  energy  and,  therefore,  would  be  expected  to  be  relatively 
more  important  for  lower  angles  of  attack.  This  is  illustrated  in  Figures  10-12,  where  the 
integrated  excitation  rate  of  the  0('D)  state  is  shown  as  a  function  of  time  for  180°,  90°  and  0°  angles 
of  attack  respectively.  The  separate  contributions  are  shown  for  the  two  mechanisms  in  the  figures, 
and  it  can  be  seen  that  the  charge  exchange  mechanism  is  calculated  to  contribute  3.2%  as  much  as 
the  collisional  excitation  mechanism  at  180°,  3.9%  at  90°,  and  12.7%  at  0°.  (These  comparisons  are 
made  at  the  end  of  the  burn.)  Hence,  the  general  trend  of  increasing  importance  of  charge 
exchange  with  decreasing  angle  of  attack  is  confirmed,  though  it  should  be  noted  that  the  charge 
exchange  rate  is  never  calculated  to  be  a  large  fraction  of  the  total  signature,  so  it  can  not  be 
confidently  deduced  that  that  portion  of  the  calculation  is  being  carried  out  properly. 

Given  the  lower  activation  energy  requirement  of  the  charge  exchange  mechanism,  it  might 
also  be  supposed  that  it  would  gain  in  importance  for  the  shutdown  transient  as  the  flow  energy 
decreases.  As  can  be  seen  in  particular  for  0°  in  Figure  12,  the  exact  opposite  Is  true  •  the  relative 
importance  of  the  charge  exchange  mechanism  drops  by  over  four  orders  of  magnitude  after 
shutdown.  The  reason  for  this  precipitous  drop  is  shown  in  Figure  13,  which  shows  the  total 
number  of  water  molecules  in  the  solution  region  as  a  function  of  time.  The  number  of  water 
molecules  understandably  drops  by  several  orders  of  magnitude  since  the  wind  and  the  exhaust 
direction  are  coaiigned.  Within  a  second  it  is  to  be  expected  that  the  vast  majority  of  the  water 
molecules,  initially  traveling  at  3.5  km/s  and  generally  accelerated  by  the  free  stream,  will  tend  to 
leave  the  solution  region  which  only  extends  four  kilometers  downstream  of  the  exit. 

The  question  remains,  however,  as  to  the  origin  of  the  collisional  excitation  aiier  the  plume 
has  largely  departed  the  solution  region.  The  explanation  for  this  is  given  in  Figure  14.  which 
shows  the  contributions  of  the  various  types  of  collisions  to  the  total  collisional  excitation  rate. 
While  the  thruster  is  on,  most  of  the  excitation  comes  from  0-N2  collisions.  This  is  to  be 
expected,  since  the  rate  constant  for  Reaction  (I)  is  quite  high,  and  N2  is  a  major  plume  species 
with  a  fairly  large  molecular  weight  (and  therefore  relative  collision  energy).  However,  aRer 
shutdown  the  curves  cross,  and  0-0  collisions  become  the  dominant  source  of  excitations.  The 
O-H2O  collisions  do  drop  in  importance  as  the  plume  leaves  the  lailution  region,  but  the  0‘N2 
collisions  do  not  drop  nearly  os  much.  This  is  because  N2  is  alst)  present  In  the  atmosphere,  so  it  is 
not  leaving  the  solution  region  in  the  same  way  that  the  H2O  is.  The  conclusion  is  that  the 
shutdown  transient  is  dominated  by  continue  collisions  within  the  atmosphere  that  remains 
disturbed  for  a  while  aRer  the  engine  shuts  down.  It  is  a  fundamental  difference  attributable  to 
looking  ti  radiation  from  an  atmospheric  rather  than  a  plume  species,  and  it  would  have  been 
difficult  to  guess  this  effect  without  a  comprehensive  model  calculation.  The  same  switch  of 
important  collision  class  can  be  seen  in  Figure  IS  for  the  180°  case,  but  the  magnitude  of  the  effect 
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is  not  as  important  here.  The  reason  is  that  for  a  retrofire  case  the  plume  molecules  remain  in  the 
solution  region  longer  since  they  are  initially  heading  into  the  free  stream  and  then  get  turned 
around  and  travel  backwards  relative  to  their  initial  motion.  The  increased  residence  time  of  plume 
species  within  the  solution  region  allows  atmosphere-plume  collisions  to  retain  their  importance  for 
a  longer  period  of  time  after  shutdown. 

The  question  of  excited  state  lifetime  is  addressed  In  Figure  16.  This  figure  shows  the 
overall  contributions  of  coliisionai  quenching  and  radiative  decay  to  the  disappearance  of  the 
0(‘D)  state  as  a  function  of  time  for  the  180°  case.  While  the  engine  is  In  operation,  an  0('D) 
atom  is  about  ten  times  as  likely  to  suffer  coliisionai  quenching  as  to  experience  radiative  decay. 
Hence,  the  effective  excited  state  lifetime  is  about  twenty  seconds  while  the  engine  is  operating,  on 
the  average.  If  just  the  brightest  portion  of  the  plume  is  inspected  (where  the  collision  and 
quenching  rates  are  particularly  high)  the  deduced  species  lifetime  is  more  like  1-2  seconds.  The 
conclusion  is  that  the  actual  species  lifetime  will  depend  greatly  on  where  in  the  plume  the  species 
is  produced.  After  shutdown,  the  two  overall  depletion  mechanisms  become  comparable,  and  the 
0(<D)  is  predicted  to  have  a  lifetime  on  the  order  of  1(X)  seconds.  Again,  this  number  would  be 
expected  to  vary  depending  on  how  coliisionai  the  region  in  which  it  is  formed  is. 
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EXCITATION  RATES  AT  180 
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Figure  10.  The  Calculated  Contributions  of  Collisional  Excitation  and  Charge 
Exchange  to  CH’D)  Production  as  a  Function  of  Tinte  for  the  1 80®  Angle  of 


EXCITATION  RATES  AT  90 
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Figure  II.  The  Calculated  Contributuins  of  Coiiisiunal  Excitation  and  Charge 
Exchange  to  0('D)  Production  as  a  Function  of  Time  for  the  90®  Angle  of 


EXCITATION  RATES  AT  0 


Figure  12.  The  Cadculaiesf  Cbnartbu!}  jas  of  CoUisional  Excitation  and  Ciiarge 
Excbai^e  to  Of’O)  Prodasiiion  as  a  Function  of  Time  for  the  0“  Angle  of 
Attack  Case. 
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Figurif  13,  The  Total  Number  of  Water  Molecules  in  the  0®  Solution  Region  as 
Funciioii  of  Time. 


COLL I SIGNAL  EXCITATIONS 
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Figure  14.  The  Contributions  of  the  Three  Types  of  Excitation  Collisions  Considered 
as  a  Function  of  Time  for  the  0°  Angle  of  Attack  Case. 


COLLISIONflL  EXCITATIONS  AT  180 


Figure  15.  The  Contributions  of  the  Three  Types  of  Excitation  Collisions  Considered 
as  a  Function  of  Time  for  the  180®  Angle  of  Attack  Case. 


REMOVAL  MECHANISMS  AT  180 
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Figure  16.  The  Competing  Effects  of  Collisional  Quenching  and  Radiative  Decay  as  a 
Function  of  Time  for  the  180®  Angle  of  Attack  Case. 


13.  SUMMARY 


The  present  version  of  the  SOCRATES  code  should  be  regarded  as  the  latest  step  in  an 
ongoing  development  project.  Particular  emphasis  has  been  paid  to  making  it  easy  for  a  non-expert 
to  use,  while  sacrificing  neither  rigor  nor  flexibility.  Future  versions  are  expected  to  allow  for  a 
gas  dynamic  interaction  with  the  solid  bodies  as  well  multiple  source  types.  Even  without  these 
features,  the  present  model  is  a  powerful  tool  which  can  be  used  for  data  analysis  and  first 
principles  predictions. 

The  sample  calculations  present  some  important  insight  into  visible  emission  from  plumes  in 
the  6300  A  region.  They  indicate  that  collisional  excitation  of  the  electronic  state  is  the  dominant 
mechanism  explaining  the  observed  emission.  This  emission  is  particularly  significant  since  it 
arises  from  the  atmosphere  itself,  so  it  is  therefore  largely  independent  of  fuel  type. 
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